Data preparation

Read in the data spreadsheets:

sheet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/LiveFishInteraction_Feb9.csv", sep = ",", strip.white = TRUE)
scored <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/ScoredFishInteraction_Feb19.csv", sep = ",", strip.white = TRUE)

# Combining the two is easier when dealing with the relevant columns alone. These are:
# Species.A, Species.B, Group.A, Group.B, and Winner.
ABwin_live <- sheet[,c(4,6,8,9,10)]
ABwin_scored <- scored[,c(10,12,13,15,16)]
both <- rbind(ABwin_live, ABwin_scored)

# Create new versions of the three data frames that include one-on-one interactions only:
sheet_oneonone <- sheet[sheet$Group.A == 1 & sheet$Group.B == 1,]
scored_oneonone <- scored[scored$Group.A == 1 & scored$Group.B == 1,]
both_oneonone <- both[both$Group.A == 1 & both$Group.B == 1,]

Clean-up

Step 1. Unidentified species are marked with an asterisk. Exclude them from the matrix based on this criterion:

drop_ufos <- function(df) {
  # Dummy vector for storing the numbers of the rows with unidentifiable species:
  tobedropped <- vector()
  for (i in 1:nrow(df)) {
    if (grepl("*", df$Species.A[i], fixed = TRUE) == TRUE |
        grepl("*", df$Species.B[i], fixed = TRUE) == TRUE |
        grepl("*", df$Winner[i], fixed = TRUE) == TRUE) {
      tobedropped <- c(tobedropped, i)
    }
  }
  return(df[-c(tobedropped),])
}

sheet2 <- drop_ufos(sheet)
scored2 <- drop_ufos(scored)
both2 <- drop_ufos(both)
sheet_ooo2 <- drop_ufos(sheet_oneonone)
scored_ooo2 <- drop_ufos(scored_oneonone)
both_ooo2 <- drop_ufos(both_oneonone)

Step 2. Remove intraspecific interactions (a few of these were scored for live observations in the beginning, before it was decided to selectively focus on interspecific interactions only):

drop_intra <- function(df) {
  tobedropped <- vector()
  for (i in 1:nrow(df)) {
    if (identical(as.character(df$Species.A[i]), as.character(df$Species.B[i])) == TRUE) {       tobedropped <- c(tobedropped, i) 
    }
  }
  if (length(tobedropped) > 0) {
    return(df[-c(tobedropped),])
  } else {
    return(df)
  }
}

sheet3 <- drop_intra(sheet2)
scored3 <- drop_intra(scored2)
both3 <- drop_intra(both2)
sheet_ooo3 <- drop_intra(sheet_ooo2)
scored_ooo3 <- drop_intra(scored_ooo2)
both_ooo3 <- drop_intra(both_ooo2)

Step 3. Drop “Striated Surgeonfish” from the matrices, as this is likely a composite category comprising both striated and brown surgeonfish. These species may have different dominance ranks, but are almost indistinguishable in appearance.

no_surgeon <- function(df) {
  return(df[df$Species.A != "Striated Surgeonfish" & df$Species.B != "Striated Surgeonfish",])
}

sheet4 <- no_surgeon(sheet3)
scored4 <- no_surgeon(scored3)
both4 <- no_surgeon(both3)
sheet_ooo4 <- no_surgeon(sheet_ooo3)
scored_ooo4 <- no_surgeon(scored_ooo3)
both_ooo4 <- no_surgeon(both_ooo3)

Generating a dominance matrix

Extract the levels of Species.A and Species.B and use them as the rows and columns of the dominance matrix:

# The union of the species listed in the "Species A" and "Species B" columns will be used
# as the rows and columns of a new matrix filled with zeros:
emptymatrix <- function(df) {
  species <- union(levels(factor(df$Species.A)), levels(factor(df$Species.B)))
  rownum <- length(species)
  newmatrix <- matrix(0, nrow = rownum, ncol = rownum)
  row.names(newmatrix) <- species
  colnames(newmatrix) <- species
  return(newmatrix)
}

dommatrix <- emptymatrix(sheet4)
scoredmatrix <- emptymatrix(scored4)
bothmatrix <- emptymatrix(both4)
domooo <- emptymatrix(sheet_ooo4)
scoredooo <- emptymatrix(scored_ooo4)
bothooo <- emptymatrix(both_ooo4)

# If species A won over species B (winner = A), add 1 to the cell in column A and row B.
# If species B won over species A (winner = B), add 1 to the cell in column B and row A.
# An increment function will be used to increase the value of each cell accordingly:
inc <- function(x) # credit: Grega Kešpret, Stack Overflow
{
 eval.parent(substitute(x <- x + 1))
}

# Now apply it to populate the matrix:
populate <- function(sheet, df) {
  for (i in 1:nrow(sheet)) {
    spec_a <- as.character(sheet$Species.A[i])
    spec_b <- as.character(sheet$Species.B[i])
    winner <- as.character(sheet$Winner[i])
    if (identical(spec_a, winner) == TRUE) {
      inc(df[as.character(spec_b), as.character(spec_a)]) 
    }
    if (identical(spec_b, winner) == TRUE) { # This should be equivalent to "else"
      inc(df[as.character(spec_a), as.character(spec_b)]) 
    }
  }
  diag(df) <- NA
  return(df)
}

dommatrix <- populate(sheet4, dommatrix)
scoredmatrix <- populate(scored4, scoredmatrix)
bothmatrix <- populate(both4, bothmatrix)
domooo <- populate(sheet_ooo4, domooo)
scoredooo <- populate(scored_ooo4, scoredooo)
bothooo <- populate(both_ooo4, bothooo)

Convert from matrices to data frames for easier manipulation:

convert_to_df <- function(df) {
  frame <- as.data.frame(df, stringsAsFactors = FALSE)
  for (i in 1:ncol(df)) {
  frame[,i] <- as.numeric(frame[,i])
  }
  return(frame)
}

domall <- convert_to_df(dommatrix)
scoredall <- convert_to_df(scoredmatrix)
bothall <- convert_to_df(bothmatrix)
liveoooall <- convert_to_df(domooo)
scoredoooall <- convert_to_df(scoredooo)
bothoooall <- convert_to_df(bothooo)

Exclude cornetfish, the only predator in the dataset. (Note: the fact that this is done here, at the matrix formatting stage rather than at the data clean-up stage, simply reflects the point at which the decision was taken during the actual analysis. Excluding a taxon from the matrix using the method below is equivalent to excluding it from the data sheet.)

# Exclude a taxon from a matrix:
# Original unvectorized function used in the pre-revisions part of the file is commented
# out below:
#
# exclude <- function(df, taxon) {
#   colrow <- c(which(row.names(df) == taxon), which(colnames(df) == taxon))
#   excluded <- df[-colrow[1], -colrow[2]]
#   return(excluded)
# }
exclude <- function(df, taxon) {
  retained <- df[!rownames(df) %in% taxon, !colnames(df) %in% taxon]
  return(retained)
}

# Note: since there are no interactions involving cornetfish in the video-only data,
# we only need to exclude it from the live-only and combined datasets.
nocornet_domall <- exclude(domall, "Cornetfish")
nocornet_bothall <- exclude(bothall, "Cornetfish")
nocornet_live <- exclude(liveoooall, "Cornetfish")
nocornet_both <- exclude(bothoooall, "Cornetfish")

Matrix completeness:

compl <- function(df) {
  counter <- 0
  for (i in 1:nrow(df)) {
    for (j in 1:ncol(df)) {
      if (!is.na(df[i,j]) & df[i,j] == 0 & df[j,i] == 0) {
        counter <- counter + 1
      }
    }
  }
  return(100*(1 - (counter/(nrow(df)*(nrow(df) - 1)))))
}
paste("The completeness of the live-only all-species matrix is ", compl(nocornet_domall), "%", sep = ""); paste("The completeness of the video-only all-species matrix is ", compl(scoredall), "%", sep = ""); paste("The completeness of the combined all-species matrix is ", compl(nocornet_bothall), "%", sep = ""); paste("The completeness of the live-only one-on-one all-species matrix is ", compl(nocornet_live), "%", sep = ""); paste("The completeness of the video-only one-on-one all-species matrix is ", compl(scoredoooall), "%", sep = ""); paste("The completeness of the combined one-on-one all-species matrix is ", compl(nocornet_both), "%", sep = "")
## [1] "The completeness of the live-only all-species matrix is 38.961038961039%"
## [1] "The completeness of the video-only all-species matrix is 27.2058823529412%"
## [1] "The completeness of the combined all-species matrix is 37.9446640316206%"
## [1] "The completeness of the live-only one-on-one all-species matrix is 33.7662337662338%"
## [1] "The completeness of the video-only one-on-one all-species matrix is 25%"
## [1] "The completeness of the combined one-on-one all-species matrix is 33.596837944664%"

Landau’s h-index:

library(gdata)
library(compete)

# To determine linearity, we need to make our dominance matrix symmetric. Below, we do
# this by subtracting the lower triangle from the upper triangle, and flipping the results
# around the diagonal:
symmatrix <- function(df) {
  # Going through the upper triangle column-wise corresponds to going through the lower 
  # triangle row-wise:
  triangle <- as.numeric(upperTriangle(df)) - as.numeric(lowerTriangle(df), byrow = TRUE)
  # Replace the upper triangle in the default column-wise order:
  symmat <- matrix(NA, nrow = nrow(df), ncol = ncol(df))
  upperTriangle(symmat) <- triangle
  # Replace the lower triangle in the row-wise order:
  lowerTriangle(symmat, byrow = TRUE) <- triangle
  return(symmat)
}

# Note that the p-values are simulation-based, and therefore stochastic: they will differ
# slightly between runs
devries(nocornet_domall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3196588 
##  p-value= 0.0015

devries(scoredall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2982113 
##  p-value= 0.0528

devries(nocornet_bothall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2962757 
##  p-value= 0.002

devries(nocornet_live, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2779799 
##  p-value= 0.0097

devries(scoredoooall, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2972647 
##  p-value= 0.0559

devries(nocornet_both, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2678731 
##  p-value= 0.0066

Uncomment the code below (and change the paths as needed) to print the matrices:

# write.table(nocornet_domall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_bothall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_live, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/live_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(scoredoooall, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/video_only_one-on-one_matrix.txt", quote = FALSE, sep = "\t")
# write.table(nocornet_both, "/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/both_datasets_one-on-one_matrix.txt", quote = FALSE, sep = "\t")

CBI

The code below was used to compute the Clutton-Brock et al. dominance index (Clutton-Brock et al. 1979, 1982).

# Remember:
# If A won over B, add 1 to [B,A] 
# If A won over C, add 1 ro [C,A]
# -> So the A column represents the species that A dominates (A wins)
# If A lost to B, add 1 to [A,B]
# If A lost to C, add 1 to [A,C]
# -> So the A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi <- function(df, species) {
  # number of entries in column 'species' that are neither 0 nor NA:
  B <- sum(df[,species] != 0, na.rm = TRUE)
  # species corresponding to those entries:
  dominates <- rownames(df)[which(df[,species] != 0)]
  # number of entries in the columns corresponding to those species that are neither 0 nor NA:
  if (length(dominates) > 0) {
    b <- sum(subset(df, select = dominates) != 0, na.rm = TRUE)
  } else {
    b <- 0
  }
  # number of entries in row B that are neither 0 nor NA:
  L <- sum(df[species,], na.rm = TRUE)
  # species corresponding to those entries:
  dominated <- colnames(df)[which(df[species,] != 0)]
  # number of entries in the rows corresponding to those species that are neither 0 nor NA:
  if (length(dominated) > 0) {
    l <- sum(subset(df, rownames(df) %in% dominated) != 0, na.rm = TRUE)
  } else {
    l <- 0
  }
  return((B + b + 1)/(L + l + 1))
}

# Apply it to the full matrix:
cbi_hierarchy <- function(df) {
  cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
  colnames(cbi_matrix) <- c("Species", "CBI")
  for (i in 1:nrow(df)) {
    cbi_matrix$Species[i] <- rownames(df)[i]
    cbi_matrix$CBI[i] <- cbi(df, rownames(df)[i])
  }
  cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
  return(cbi_matrix)
}

cbi_ncdomall <- cbi_hierarchy(nocornet_domall)
cbi_scoredall <- cbi_hierarchy(scoredall)
cbi_ncbothall <- cbi_hierarchy(nocornet_bothall)
cbi_ncooolive <- cbi_hierarchy(nocornet_live)
cbi_ooovideo <- cbi_hierarchy(scoredoooall)
cbi_ncoooboth <- cbi_hierarchy(nocornet_both)

Frequency of winning

win_freq <- function(df) {
  j <- 1
  species <- union(levels(df$Species.A), levels(df$Species.B))
  freq_table <- data.frame(matrix(NA, nrow = length(species), ncol = 2))
  colnames(freq_table) <- c("Species", "WinFreq")
  for(i in species) {
    if (grepl("*", i, fixed = TRUE) == TRUE) {
      next
    } else {
      freq_table[j,1] <- i
      involved <- sum(df$Species.A == i) + sum(df$Species.B == i)
      won <- sum(df$Winner == i)
      freq_table[j,2] <- won/involved
      j <- j + 1
    }
  }
  return(freq_table[complete.cases(freq_table),])
}

Sizes

Compare the sizes estimated from live observations with those reported by Randall (2005):

# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", strip.white = TRUE, header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]

# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)

# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
  return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
  return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])

# Do book sizes correlate with live observation sizes?
# a <- cor.test(meansizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(a$estimate[[1]],5), ", p-value: ", a$p.value, sep = "")
# b <- cor.test(meansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(b$estimate[[1]],5), ", p-value: ", b$p.value, sep = "")
# c <- cor.test(mediansizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(c$estimate[[1]],5), ", p-value: ", c$p.value, sep = "")
d <- cor.test(meansizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(d$estimate[[1]],5), ", p-value: ", d$p.value, sep = "")
## [1] "correlation: 0.9969, p-value: 1.4440712484574e-23"
e <- cor.test(meansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(e$estimate[[1]],5), ", p-value: ", e$p.value, sep = "")
## [1] "correlation: 0.94797, p-value: 2.11724427906556e-11"
f <- cor.test(mediansizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(f$estimate[[1]],5), ", p-value: ", f$p.value, sep = "")
## [1] "correlation: 0.94834, p-value: 1.973173161907e-11"
# Non-unique ranks, gives an error:
# g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")
# h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")
# i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")
# paste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")
j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")
## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"
k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")
## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"
l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")
## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"
# Reassign the mediansizes vector so that it includes Juvenile Parrotfish:
mediansizes <- colMedians(fishsize[,c(2:4,6:25)])

# Create a new data frame that associates each species with the median of its observed
# size distribution, excluding Cornetfish and "Striated Surgeonfish":
obs_sizes <- data.frame(matrix(ncol = 2, nrow = length(mediansizes)))
colnames(obs_sizes) <- c("Species", "Size")
obs_sizes[,1] <- as.character(sizes$Species[c(1:3,5:24)])
obs_sizes[,2] <- as.numeric(mediansizes)

Include mass (scalled isometrically and allometrically) as an additional variable:

mass <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/body_to_mass_data.csv", sep = ",", strip.white = TRUE)
colnames(mass)[1] <- "Species"
lengthmass <- merge(obs_sizes, mass, by = "Species")
lengthmass[,5] <- lengthmass[,3]*(lengthmass[,2]^lengthmass[,4])
lengthmass[,6] <- lengthmass[,5]^(2/3)
colnames(lengthmass)[c(2,5,6)] <- c("Length", "Mass", "Mass2_3")

Test of mass vs. individual fighting ability

Test the relationship between the two measures of dominance (CBI and winning frequency) and body mass or body mass raised to 2/3 by creating a single data frame that contains all the four variables. Note: from this point on, we focus exclusively on the combined (live plus video) dataset of 23 species (no cornetfish) and one-on-one-interactions.

both_test <- merge(cbi_ncoooboth, win_freq(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",]), by = "Species")
both_test <- merge(both_test, lengthmass, by = "Species")
# Number of interactions:
nrow(both_ooo4[both_ooo4$Species.A != "Cornetfish" & both_ooo4$Species.B != "Cornetfish",])
## [1] 401
cbi_mass <- lm(both_test$CBI ~ both_test$Mass)
summary(cbi_mass)$coefficients[2, 4]
## [1] 0.4754628
summary(cbi_mass)$r.squared
## [1] 0.02452714
winfreq_mass <- lm(both_test$WinFreq ~ both_test$Mass)
summary(winfreq_mass)$coefficients[2, 4]
## [1] 0.1346416
summary(winfreq_mass)$r.squared
## [1] 0.1033787
cbi_mass2_3 <- lm(both_test$CBI ~ both_test$Mass2_3)
summary(cbi_mass2_3)$coefficients[2, 4]
## [1] 0.2865715
summary(cbi_mass2_3)$r.squared
## [1] 0.05387003
winfreq_mass2_3 <- lm(both_test$WinFreq ~ both_test$Mass2_3)
summary(winfreq_mass2_3)$coefficients[2, 4]
## [1] 0.1029628
summary(winfreq_mass2_3)$r.squared
## [1] 0.1215855
# How far is spotted porcupinefish from the median in terms of the interquartile range?
# (Note: this is a simple measure of the extremeness of an outlier)
(both_test$Mass[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass))/IQR(both_test$Mass)
## [1] 9.456649
(both_test$Mass2_3[both_test$Species == "Spotted Porcupinefish"] - median(both_test$Mass2_3))/IQR(both_test$Mass2_3)
## [1] 4.972994
# Exclude spotted porcupinefish:
nocornet_nosp_both <- exclude(nocornet_both, "Spotted Porcupinefish")
cbi_ncoooboth_nosp <- cbi_hierarchy(nocornet_nosp_both)
both_test_nosp <- merge(cbi_ncoooboth_nosp, win_freq(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),]), by = "Species")
both_test_nosp <- merge(both_test_nosp, lengthmass, by = "Species")
# Number of data points:
nrow(both_ooo4[!both_ooo4$Species.A %in% c("Cornetfish", "Spotted Porcupinefish") & !both_ooo4$Species.B %in% c("Cornetfish", "Spotted Porcupinefish"),])
## [1] 398
cbi_mass_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass)
summary(cbi_mass_nosp)$coefficients[2, 4]
## [1] 0.0301074
summary(cbi_mass_nosp)$r.squared
## [1] 0.214146
winfreq_mass_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass)
summary(winfreq_mass_nosp)$coefficients[2, 4]
## [1] 0.2470032
summary(winfreq_mass_nosp)$r.squared
## [1] 0.06638793
cbi_mass2_3_nosp <- lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3)
summary(cbi_mass2_3_nosp)$coefficients[2, 4]
## [1] 0.03673997
summary(cbi_mass2_3_nosp)$r.squared
## [1] 0.2003125
winfreq_mass2_3_nosp <- lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3)
summary(winfreq_mass2_3_nosp)$coefficients[2, 4]
## [1] 0.2267295
summary(winfreq_mass2_3_nosp)$r.squared
## [1] 0.07216368

This is the code used to generate Figure 1 of the manuscript (NOTE: This figure summarizes the results of the regression analyses uncorrected for phylogenetic relatedness, and is now given as Figure S3 in the Supplementary Information.)

# png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Supp_Figure_3.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))

par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(lm(both_test$CBI ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass)$r.squared
pvalue <- summary(cbi_mass)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass_nosp)$coefficients[2,4]
text(800, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$CBI ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$CBI ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(cbi_mass2_3)$r.squared
pvalue <- summary(cbi_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(cbi_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(cbi_mass2_3_nosp)$coefficients[2,4]
text(15, 10.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(15, 9.1, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(15, 7.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(15, 6, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_test$Mass, both_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass)$r.squared
pvalue <- summary(winfreq_mass)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass_nosp)$coefficients[2,4]
text(800, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(800, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(800, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(800, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_test$Mass2_3, both_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_test$Mass == both_test$Mass[both_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(lm(both_test$WinFreq ~ both_test$Mass2_3), lwd = 5)
abline(lm(both_test_nosp$WinFreq ~ both_test_nosp$Mass2_3), lwd = 5, col = "gray70")
rsquared <- summary(winfreq_mass2_3)$r.squared
pvalue <- summary(winfreq_mass2_3)$coefficients[2,4]
rsquared_nosp <- summary(winfreq_mass2_3_nosp)$r.squared
pvalue_nosp <- summary(winfreq_mass2_3_nosp)$coefficients[2,4]
text(60, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared, 3))), pos = 4)
text(60, 0.4 - 1.3/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(60, 0.4 - (1.3+1.8)/16, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(rsquared_nosp, 3))), col = "gray70", pos = 4)
text(60, 0.4 - (2*1.3+1.8)/16, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()

Testing Lanchester’s Laws

For each interaction, multiply the body mass of either species (scaled isometrically or allometrically) by its group size or group size squared. See which species maximizes the product, and whether it matches the observed winner. Run a binomial test on the proportion of successful predictions.

both_group <- both4[both4$Group.A != both4$Group.B,]

# Clean-up: drop rows containing NAs:
both_group <- both_group[complete.cases(both_group),]

# Exclude cornetfish:
both_group <- both_group[both_group$Species.A != "Cornetfish" &
                         both_group$Species.B != "Cornetfish",]

# Number of data points:
nrow(both_group)
## [1] 83
# Get the body sizes for both competitors:
for(i in 1:nrow(both_group)) {
  species_a <- as.character(both_group$Species.A[i])
  species_b <- as.character(both_group$Species.B[i])
  both_group[i,6] <- lengthmass$Mass[lengthmass$Species == species_a]
  both_group[i,7] <- lengthmass$Mass[lengthmass$Species == species_b]
}
colnames(both_group)[6:7] <- c("BodySize.A", "BodySize.B")

# Convert group sizes to the right object type (from 'character' to 'double'):
both_group$Group.A <- as.numeric(both_group$Group.A)
both_group$Group.B <- as.numeric(both_group$Group.B)

# Make predictions based on the linear law:
for(i in 1:nrow(both_group)) {
  if (both_group$Group.A[i]*both_group$BodySize.A[i] > both_group$Group.B[i]*both_group$BodySize.B[i]) {
    both_group[i,8] <- both_group$Species.A[i] 
  } else {
    both_group[i,8] <- both_group$Species.B[i] 
  }
}

# Make predictions based on the square law:
for(i in 1:nrow(both_group)) {
  if ((both_group$Group.A[i]^2)*both_group$BodySize.A[i] > (both_group$Group.B[i]^2)*both_group$BodySize.B[i]) {
    both_group[i,9] <- both_group$Species.A[i] 
  } else {
    both_group[i,9] <- both_group$Species.B[i] 
  }
}
colnames(both_group)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")

# Test if the linear law works:
for(i in 1:nrow(both_group)) {
  if (identical(as.character(both_group$Winner[i]), as.character(both_group$Linear.Law.Prediction[i])) == TRUE) {
    both_group[i,10] <- 1
  } else {
    both_group[i,10] <- 0
  }
}

# Test if the square law works:
for(i in 1:nrow(both_group)) {
  if (identical(as.character(both_group$Winner[i]), as.character(both_group$Square.Law.Prediction[i])) == TRUE) {
    both_group[i,11] <- 1
  } else {
    both_group[i,11] <- 0
  }
}
colnames(both_group)[10:11] <- c("Linear.True", "Square.True")

# Proportion of successful predictions:
sum(both_group$Linear.True)/nrow(both_group)
## [1] 0.6626506
sum(both_group$Square.True)/nrow(both_group)
## [1] 0.4939759
# Test significance:
dbinom(sum(both_group$Linear.True), size = nrow(both_group), prob = 1/2) 
## [1] 0.001053886
dbinom(sum(both_group$Square.True), size = nrow(both_group), prob = 1/2)
## [1] 0.08679764
# When do the two laws give different predictions?
dif <- both_group[both_group$Linear.Law.Prediction != both_group$Square.Law.Prediction,]
dbinom(sum(dif$Linear.True), size = nrow(dif), prob = 1/2) 
## [1] 0.001087189
# Now repeat all of this for body mass raised to two thirds:
both_group2 <- both4[both4$Group.A != both4$Group.B,]
both_group2 <- both_group2[complete.cases(both_group2),]
both_group2 <- both_group2[both_group2$Species.A != "Cornetfish" &
                           both_group2$Species.B != "Cornetfish",]

for(i in 1:nrow(both_group2)) {
  species_a <- as.character(both_group2$Species.A[i])
  species_b <- as.character(both_group2$Species.B[i])
  both_group2[i,6] <- lengthmass$Mass2_3[lengthmass$Species == species_a]
  both_group2[i,7] <- lengthmass$Mass2_3[lengthmass$Species == species_b]
}
colnames(both_group2)[6:7] <- c("BodySize.A", "BodySize.B")

both_group2$Group.A <- as.numeric(both_group2$Group.A)
both_group2$Group.B <- as.numeric(both_group2$Group.B)

for(i in 1:nrow(both_group2)) {
  if (both_group2$Group.A[i]*both_group2$BodySize.A[i] > both_group2$Group.B[i]*both_group2$BodySize.B[i]) {
    both_group2[i,8] <- both_group2$Species.A[i] 
  } else {
    both_group2[i,8] <- both_group2$Species.B[i] 
  }
}

for(i in 1:nrow(both_group2)) {
  if ((both_group2$Group.A[i]^2)*both_group2$BodySize.A[i] > (both_group2$Group.B[i]^2)*both_group2$BodySize.B[i]) {
    both_group2[i,9] <- both_group2$Species.A[i] 
  } else {
    both_group2[i,9] <- both_group2$Species.B[i] 
  }
}
colnames(both_group2)[8:9] <- c("Linear.Law.Prediction", "Square.Law.Prediction")

for(i in 1:nrow(both_group2)) {
  if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Linear.Law.Prediction[i])) == TRUE) {
    both_group2[i,10] <- 1
  } else {
    both_group2[i,10] <- 0
  }
}

for(i in 1:nrow(both_group2)) {
  if (identical(as.character(both_group2$Winner[i]), as.character(both_group2$Square.Law.Prediction[i])) == TRUE) {
    both_group2[i,11] <- 1
  } else {
    both_group2[i,11] <- 0
  }
}
colnames(both_group2)[10:11] <- c("Linear.True", "Square.True")

sum(both_group2$Linear.True)/nrow(both_group2)
## [1] 0.5662651
sum(both_group2$Square.True)/nrow(both_group2)
## [1] 0.3614458
dbinom(sum(both_group2$Linear.True), size = nrow(both_group2), prob = 1/2) 
## [1] 0.04240454
dbinom(sum(both_group2$Square.True), size = nrow(both_group2), prob = 1/2)
## [1] 0.003597748
dif2 <- both_group2[both_group2$Linear.Law.Prediction != both_group2$Square.Law.Prediction,]
dbinom(sum(dif2$Linear.True), size = nrow(dif2), prob = 1/2) 
## [1] 0.0001001358

Supplementary analyses

“Rarefaction” analysis

Plot a curve showing how many unique interactions you get with each new video scored:

curve_analysis <- drop_ufos(scored)
curve_analysis <- curve_analysis[complete.cases(curve_analysis),]
curve_analysis <- curve_analysis[,c(4,5,6,10,13)]
# Convert the dates and times into a readable format and sort the dataset by them:
library(chron)
curve_analysis$Date.of.recording <- as.Date(curve_analysis$Date.of.recording, format = "%d %B %y")
# This will makes the times readable by 'chron':
curve_analysis$Time.of.Recording <- paste(curve_analysis$Time.of.Recording, ":00", sep = "")
curve_analysis$Time.of.Recording <- chron(times = curve_analysis$Time.of.Recording)
# Sort by date and time:
curve_analysis <- curve_analysis[order(curve_analysis$Date.of.recording, curve_analysis$Time.of.Recording, decreasing = TRUE),]

species_pairs <- list()
rarefaction <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis)) {
  new_pair <- c(as.character(curve_analysis$Species.A[i]), as.character(curve_analysis$Species.B[i]))
  species_pairs <- c(species_pairs, list(new_pair))
  rarefaction[i,1] <- i
  rarefaction[i,2] <- length(unique(species_pairs))
}
plot(rarefaction$Observation, rarefaction$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species")

# Now order the observations according to scoring time rather than recording time:
curve_analysis2 <- drop_ufos(scored)
curve_analysis2 <- curve_analysis2[complete.cases(curve_analysis2),]
curve_analysis2 <- curve_analysis2[,c(3,10,13)]
curve_analysis2$Date.of.scoring <- as.Date(curve_analysis2$Date.of.scoring, format = "%d %B %Y")
curve_analysis2 <- curve_analysis2[order(curve_analysis2$Date.of.scoring, decreasing = TRUE),]

species_pairs2 <- list()
rarefaction2 <- data.frame(matrix(NA, nrow = nrow(curve_analysis), ncol = 2))
colnames(rarefaction2) <- c("Observation", "Unique Pairs")
for(i in 1:nrow(curve_analysis2)) {
  new_pair <- c(as.character(curve_analysis2$Species.A[i]), as.character(curve_analysis2$Species.B[i]))
  species_pairs2 <- c(species_pairs2, list(new_pair))
  rarefaction2[i,1] <- i
  rarefaction2[i,2] <- length(unique(species_pairs2))
}
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")

This is the code used to generate Figure S4:

png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.png", width = 2880, height = 1800, pointsize = 64)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_S1.eps")
plot(rarefaction2$Observation, rarefaction2$`Unique Pairs`, pch = 16, xlab = "Observations", ylab = "Unique pairs of interacting species", main = "Contribution of new observations to dataset expansion")
dev.off()

Breakdown by diet

diet <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/diet.csv", sep = ",", strip.white = TRUE)

# Add dietary information for both Species A and Species B
for(i in 1:nrow(both_group)) {
  both_group[i, 12] <- diet$Diet[diet$Species == as.character(both_group$Species.A[i])]
  both_group[i, 13] <- diet$Diet[diet$Species == as.character(both_group$Species.B[i])]
}
colnames(both_group)[12:13] <- c("Diet.A", "Diet.B")

# Now subsample the test frame to include only those interactions that take place between
# species with the same diet:
test_benth <- both_group[both_group$Diet.A == "Benthic Invertebrate Consumer" &
                         both_group$Diet.B == "Benthic Invertebrate Consumer",]
nrow(test_benth)
## [1] 0
test_omni <- both_group[both_group$Diet.A == "Omnivore" & both_group$Diet.B == "Omnivore",]
nrow(test_omni)
## [1] 16
test_herb <- both_group[both_group$Diet.A == "Herbivore/Detritivore" & both_group$Diet.B == "Herbivore/Detritivore",]
nrow(test_herb)
## [1] 4
test_coral <- both_group[both_group$Diet.A == "Corallivore" & both_group$Diet.B == "Corallivore",]
nrow(test_coral)
## [1] 0
# Test the validity of Lanchester's laws for interaction subsets of non-zero length.
# Omnivores:
sum(test_omni$Linear.True)/nrow(test_omni)
## [1] 0.8125
sum(test_omni$Square.True)/nrow(test_omni)
## [1] 0.5625
dbinom(sum(test_omni$Linear.True), size = nrow(test_omni), prob = 1/2) 
## [1] 0.008544922
dbinom(sum(test_omni$Square.True), size = nrow(test_omni), prob = 1/2)
## [1] 0.1745605
dif3 <- test_omni[test_omni$Linear.Law.Prediction != test_omni$Square.Law.Prediction,]
nrow(dif3)
## [1] 4
sum(dif3$Linear.True)/nrow(dif3)
## [1] 1
dbinom(sum(dif3$Linear.True), size = nrow(dif3), prob = 1/2) 
## [1] 0.0625
# Herbivores/Detritivores:
sum(test_herb$Linear.True)/nrow(test_herb)
## [1] 1
sum(test_herb$Square.True)/nrow(test_herb)
## [1] 0.25
dbinom(sum(test_herb$Linear.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625
dbinom(sum(test_herb$Square.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.25
dif4 <- test_herb[test_herb$Linear.Law.Prediction != test_herb$Square.Law.Prediction,]
nrow(dif4)
## [1] 3
sum(dif4$Linear.True)/nrow(dif4)
## [1] 1
dbinom(sum(dif4$Linear.True), size = nrow(dif4), prob = 1/2) 
## [1] 0.125
# Re-run everything with mass raised to two thirds:
for(i in 1:nrow(both_group2)) {
  both_group2[i, 12] <- diet$Diet[diet$Species == as.character(both_group2$Species.A[i])]
  both_group2[i, 13] <- diet$Diet[diet$Species == as.character(both_group2$Species.B[i])]
}
colnames(both_group2)[12:13] <- c("Diet.A", "Diet.B")
test_benth2 <- both_group2[both_group2$Diet.A == "Benthic Invertebrate Consumer" &
                         both_group2$Diet.B == "Benthic Invertebrate Consumer",]
test_omni2 <- both_group2[both_group2$Diet.A == "Omnivore" & both_group2$Diet.B == "Omnivore",]
test_herb2 <- both_group2[both_group2$Diet.A == "Herbivore/Detritivore" & both_group2$Diet.B == "Herbivore/Detritivore",]

sum(test_omni2$Linear.True)/nrow(test_omni2)
## [1] 0.5625
sum(test_omni2$Square.True)/nrow(test_omni2)
## [1] 0.25
dbinom(sum(test_omni2$Linear.True), size = nrow(test_omni2), prob = 1/2) 
## [1] 0.1745605
dbinom(sum(test_omni2$Square.True), size = nrow(test_omni2), prob = 1/2) 
## [1] 0.027771
dif5 <- test_omni2[test_omni2$Linear.Law.Prediction != test_omni2$Square.Law.Prediction,]
nrow(dif5)
## [1] 5
sum(dif5$Linear.True)/nrow(dif5)
## [1] 1
dbinom(sum(dif5$Linear.True), size = nrow(dif5), prob = 1/2) 
## [1] 0.03125
sum(test_herb2$Linear.True)/nrow(test_herb2)
## [1] 0.75
sum(test_herb2$Square.True)/nrow(test_herb2)
## [1] 0
dbinom(sum(test_herb2$Linear.True), size = nrow(test_herb2), prob = 1/2)
## [1] 0.25
dbinom(sum(test_herb2$Square.True), size = nrow(test_herb2), prob = 1/2)
## [1] 0.0625
dif6 <- test_herb2[test_herb2$Linear.Law.Prediction != test_herb2$Square.Law.Prediction,]
nrow(dif6)
## [1] 3
sum(dif6$Linear.True)/nrow(dif6)
## [1] 1
dbinom(sum(dif6$Linear.True), size = nrow(dif6), prob = 1/2) 
## [1] 0.125

Analyses not included in the manuscript

These include exploratory analyses as well as applications of previously proposed methods to the fish dataset (e.g., Shelley et al.’s (2004) logistic regression).

Most often interacting species

Subsample the data sheets down to the 10 most often interacting species. This will increase completeness (substantially) and linearity (slightly), but give you fewer data points to work with.

mostoften <- function(sheet, df, n) {
  occurrences <- tail(sort(table(unlist(sheet[,c("Species.A", "Species.B")]))), n)
  max_n <- as.character(as.data.frame(occurrences)[,1])
  max_n_matrix <- subset(df, select = max_n)
  max_n_matrix <- subset(max_n_matrix, rownames(max_n_matrix) %in% max_n)
  max_n_matrix <- max_n_matrix[order(row.names(max_n_matrix)), order(colnames(max_n_matrix))]
  return(max_n_matrix)
}

dommatrix_max10 <- mostoften(sheet4, dommatrix, 10)
scoredmatrix_max10 <- mostoften(scored4, scoredmatrix, 10)
bothmatrix_max10 <- mostoften(both4, bothmatrix, 10)
domooo_max10 <- mostoften(sheet_ooo4, domooo, 10)
scoredooo_max10 <- mostoften(scored_ooo4, scoredooo, 10)
bothooo_max10 <- mostoften(both_ooo4, bothooo, 10)

dommax10 <- convert_to_df(dommatrix_max10)
scoredmax10 <- convert_to_df(scoredmatrix_max10)
bothmax10 <- convert_to_df(bothmatrix_max10)
liveooomax10 <- convert_to_df(domooo_max10)
scoredooomax10 <- convert_to_df(scoredooo_max10)
bothooomax10 <- convert_to_df(bothooo_max10)

paste("The completeness of the live-only 10-species matrix is ", compl(dommax10), "%", sep = "")
paste("The completeness of the video-only 10-species matrix is ", compl(scoredmax10), "%", sep = "")
paste("The completeness of the combined 10-species matrix is ", compl(bothmax10), "%", sep = "")
paste("The completeness of the live-only one-on-one 10-species matrix is ", compl(liveooomax10), "%", sep = "")
paste("The completeness of the video-only one-on-one 10-species matrix is ", compl(scoredooomax10), "%", sep = "")
paste("The completeness of the combined one-on-one 10-species matrix is ", compl(bothooomax10), "%", sep = "")

cbi_max10 <- cbi_hierarchy(dommax10)
cbi_scoredmax10 <- cbi_hierarchy(scoredmax10)
cbi_bothmax10 <- cbi_hierarchy(bothmax10)
cbi_liveooomax10 <- cbi_hierarchy(liveooomax10)
cbi_scoredooomax10 <- cbi_hierarchy(scoredooomax10)
cbi_bothooomax10 <- cbi_hierarchy(bothooomax10)

devries(dommax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothmax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(liveooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(scoredooomax10, Nperms = 10000, history = FALSE, plot = TRUE)
devries(bothooomax10, Nperms = 10000, history = FALSE, plot = TRUE)

I&SI rank

A third measure of dominance that is particularly applicable to weakly but significantly linear dominance hierarchies, implemented as described in Schmid & de Vries (2013).

library(compete)
i_and_si_liveooo <- isi13(liveoooall, nTries = 10)
rev(i_and_si_liveooo$best_order) # The order is given from least to most dominant by default

i_and_si_videoooo <- isi13(scoredoooall, nTries = 10)
rev(i_and_si_videoooo$best_order)

i_and_si_bothooo <- isi13(bothoooall, nTries = 10)
rev(i_and_si_bothooo$best_order)

i_and_si <- isi13(dommax10, nTries = 10)
rev(i_and_si$best_order) 

i_and_si_scored <- isi13(scoredmax10, nTries = 10)
rev(i_and_si_scored$best_order)

i_and_si_both <- isi13(bothmax10, nTries = 10)
rev(i_and_si_both$best_order)

# Create a data frame with alphabetically ordered species and their corresponding I&SI ranks:
isi_both <- i_and_si_bothooo$best_order
isi_both <- setdiff(isi_both, "Cornetfish")
test_isi <- data.frame(matrix(NA, nrow = length(isi_both), ncol = 2))
colnames(test_isi) <- c("Species", "I_SI_Rank")
test_isi$Species <- sort(isi_both)
test_isi$I_SI_Rank <- match(sort(isi_both), isi_both)

Nonlinear models

Do simple nonlinear models (exponential, logarithmic, quadratic) explain the relationship between individual fighting ability (as determined by CBI and the frequency of winning) and body mass?

# Add mass squared to the test frame:
both_test[,9] <- both_test$Mass^2
colnames(both_test)[9] <- "Mass_Squared"

stats <- function(df) {
  results <- data.frame(matrix(NA, nrow = 2, ncol = 14))
  # Columns represent different statistics associated with each model. Spearman's rho
  # and Pearson's coefficient are included for comparison; if the former is significant
  # but the latter is not, this shows that the relationship is monotonic but nonlinear.
  colnames(results) <- c("rho", "Spear_p-value", "cor", "Pear_p-value", "linAIC", "expRsq", "exp_p-value", "expAIC", "logRsq", "log_p-value", "logAIC", "quadrRsq", "quadr_p-value", "quadrAIC")
  # Rows represent dependent variables:
  row.names(results) <- c("CBI_Mass", "WinFreq_Mass")
  spear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "spearman", alternative = "two.sided")
  spear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "spearman", alternative = "two.sided")
  pear_mass_CBI <- cor.test(df$CBI, df$Mass, method = "pearson", alternative = "two.sided")
  pear_mass_WinFreq <- cor.test(df$WinFreq, df$Mass, method = "pearson", alternative = "two.sided")
  linear_mass_CBI <- lm(df$CBI ~ df$Mass)
  expo_mass_CBI <- lm(log(df$CBI) ~ df$Mass)
  loga_mass_CBI <- lm(df$CBI ~ log(df$Mass))
  quadra_mass_CBI <- lm(df$CBI ~ df$Mass + df$Mass_Squared)
  linear_mass_WinFreq <- lm(df$WinFreq  ~ df$Mass)
  # No exponential model for winning frequency: the minimum winning frequency observed
  # in the dataset is 0, and simple log transformation therefore cannot be applied to it.
  loga_mass_WinFreq <- lm(df$WinFreq  ~ log(df$Mass))
  quadra_mass_WinFreq <- lm(df$WinFreq  ~ df$Mass + df$Mass_Squared)
  results[1,1] <- spear_mass_CBI$estimate; results[1,2] <- spear_mass_CBI$p.value
  results[1,3] <- pear_mass_CBI$estimate; results[1,4] <- pear_mass_CBI$p.value
  results[1,5] <- AIC(linear_mass_CBI); results[1,6] <- summary(expo_mass_CBI)$r.squared
  results[1,7] <- summary(expo_mass_CBI)$coefficients[,4][[2]]
  results[1,8] <- AIC(expo_mass_CBI); results[1,9] <- summary(loga_mass_CBI)$r.squared
  results[1,10] <- summary(loga_mass_CBI)$coefficients[,4][[2]]
  results[1,11] <- AIC(loga_mass_CBI); results[1,12] <- summary(quadra_mass_CBI)$r.squared
  results[1,13] <- summary(quadra_mass_CBI)$coefficients[,4][[2]]
  results[1,14] <- AIC(quadra_mass_CBI)
  results[2,1] <- spear_mass_WinFreq$estimate; results[2,2] <- spear_mass_WinFreq$p.value
  results[2,3] <- pear_mass_WinFreq$estimate; results[2,4] <- pear_mass_WinFreq$p.value
  results[2,5] <- AIC(linear_mass_WinFreq)
  results[2,9] <- summary(loga_mass_WinFreq)$r.squared
  results[2,10] <- summary(loga_mass_WinFreq)$coefficients[,4][[2]]
  results[2,11] <- AIC(loga_mass_WinFreq)
  results[2,12] <- summary(quadra_mass_WinFreq)$r.squared
  results[2,13] <- summary(quadra_mass_WinFreq)$coefficients[,4][[2]]
  results[2,14] <- AIC(quadra_mass_WinFreq)
  return(results)
}

both_test_results <- stats(both_test)

# Adjust the p-values for multiple comparisons:
library(stats)
pvalues <- unlist(both_test_results[,c(2,4,7,10,13)])
p.adjust(pvalues, "bonferroni", sum(!is.na(pvalues)))

Logistic regression

The code below is an implementation of the logistic regression analysis described by Shelley et al. (2004). Given the paucity of data, no pairwise comparisons were made; each of the selected species was compared against all other species it interacted with. Unlike Shelley et al. (2004), we regress against both group size and group size squared.

# Extract those species that participated in an interaction with more than one individual
# at least a given number of times:
atleastoncemorethan <- function(df, n) {
  output <- vector()
  species <- union(unique(df$Species.A), unique(df$Species.B))
  for (i in species) {
    a.frame <- df[df$Species.A == i,]
    b.frame <- df[df$Species.B == i,]
    if (sum(!is.na(a.frame$Group.A) & a.frame$Group.A > 1) + sum(!is.na(b.frame$Group.B) & b.frame$Group.B > 1) >= n) {
      output <- c(output, as.character(i))
    }
  }
  return(unique(output))
}

# Some ad hoc data clean-up before applying the atleastoncemorethan() function:
scored4$Group.B <- as.numeric(as.character(scored4$Group.B))
both4$Group.B <- as.numeric(both4$Group.B)

# Include only those species that had a group size greater than 1 at least twice. This
# ensures that there is variability in the predictor variable.
live_logres <- sheet4[sheet4$Species.A %in% atleastoncemorethan(sheet4, 2) &
                      sheet4$Species.B %in% atleastoncemorethan(sheet4, 2),]
video_logres <- scored4[scored4$Species.A %in% atleastoncemorethan(scored4, 2) &
                        scored4$Species.B %in% atleastoncemorethan(scored4, 2),]
both_logres <- both4[both4$Species.A %in% atleastoncemorethan(both4, 2) &
                    both4$Species.B %in% atleastoncemorethan(both4, 2),]

# Convert from wide to long format: every two successive rows correspond to one
# interaction, with the winner coded as "1" and the loser coded as "0".
wide_to_long <- function(df) {
  longframe <- data.frame(matrix(NA, ncol = 3, nrow = 2*nrow(df)))
  colnames(longframe) <- c("Species", "Group_Size", "Is_Winner")
  for (i in 1:nrow(df)) {
    longframe$Species[2*i-1] <- as.character(df$Species.A[i])
    longframe$Species[2*i] <- as.character(df$Species.B[i])
    longframe$Group_Size[2*i-1] <- df$Group.A[i]
    longframe$Group_Size[2*i] <- df$Group.B[i]
    if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
      longframe$Is_Winner[2*i-1] <- 1
      longframe$Is_Winner[2*i] <- 0
    }
    if (identical(as.character(df$Species.B[i]), as.character(df$Winner[i])) == TRUE) {
      longframe$Is_Winner[2*i-1] <- 0
      longframe$Is_Winner[2*i] <- 1
    }
  }
  return(longframe)
}

# Make sure winner is always either A or B, and that you have group sizes for all rows:
live_long <- wide_to_long(live_logres)
nrow(live_long) == nrow(live_long[!is.na(live_long$Is_Winner),]) 
nrow(live_long) == nrow(live_long[!is.na(live_long$Group_Size),])
video_long <- wide_to_long(video_logres)
nrow(video_long) == nrow(video_long[!is.na(video_long$Is_Winner),])
nrow(video_long) == nrow(video_long[!is.na(video_long$Group_Size),])
both_long <- wide_to_long(both_logres)
nrow(both_long) == nrow(both_long[!is.na(both_long$Is_Winner),])
nrow(both_long) == nrow(both_long[!is.na(both_long$Group_Size),])

# Add group size squared
square_it <- function(df) {
  a <- ncol(df)
  df[,a+1] <- df$Group_Size^2
  colnames(df)[a+1] <- "Group_Size_Squared"
  return(df)
}

live_long <- square_it(live_long)
video_long <- square_it(video_long)
both_long <- square_it(both_long)

# Extract those species that won and lost at least 1 interaction. This ensures that there
# is variability in the dependent variable:
won_once <- function(df) {
  output <- vector()
  species <- unique(df$Species)
  for(i in species) {
    subset.frame <- df[df$Species == i,]
    if(0 %in% subset.frame$Is_Winner & 1 %in% subset.frame$Is_Winner) {
      output <- c(output, i)
    }
  }
  return(df[df$Species %in% output,])
}

# This function creates a data frame with logistic regression results. The regression is
# run against group size if x is FALSE and against group size squared if x is TRUE.
library(pscl)
species_logistic <- function(df, x = TRUE) {
  j <- 1
  species <- unique(df$Species)
  results <- data.frame(matrix(NA, nrow = length(species), ncol = 6))
  colnames(results) <- c("Species", "Coefficient", "p-value", "McFaddenpseudoRsq", "AIC", "n")
  for(i in species) {
    if (x == FALSE) {
      model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size[df$Species == i])
    } else {
      model <- glm(df$Is_Winner[df$Species == i] ~ df$Group_Size_Squared[df$Species == i] + df$Group_Size[df$Species == i])
    }
    results[j,1] <- i
    results[j,2] <- model$coefficients[2]
    results[j,3] <- coef(summary(model))[,4][2]
    results[j,4] <- pR2(model)[4]
    results[j,5] <- AIC(model)
    results[j,6] <- nrow(df[df$Species == i,])
    j <- j + 1
  }
  return(results)
}

# Linear group size
species_logistic(won_once(live_long), FALSE)
species_logistic(won_once(video_long), FALSE)
species_logistic(won_once(both_long), FALSE)

# Group size squared
species_logistic(won_once(live_long), TRUE)
species_logistic(won_once(video_long), TRUE)
species_logistic(won_once(both_long), TRUE)

Chi-squared

This is an alternative way of testing Lanchester’s laws. Instead of applying the binomial test to a binary variable (prediction: correct, incorrect), we apply \(\chi^2\) to a 4-entry contingency table (correctly predicted win, incorrectly predicted win, correctly predicted loss, incorrectly predicted loss):

populate_cont_table <- function(df, law) {
  # Create a contingency table whose rows represent (1) Wins, (2) Losses predicted
  # by the linear law and whose columns represent observed (1) Wins and (2) Losses
  cont_table <- matrix(0, nrow = 2, ncol = 2)
  # Focus solely on Species A. This is arbitrary; we could instead focus on B by
  # interchanging individual rows and individual columns.
  for(i in 1:nrow(df)) {
    # A won:
    if (identical(as.character(df$Species.A[i]), as.character(df$Winner[i])) == TRUE) {
      # A was predicted to win:
      if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
        inc(cont_table[1,1])
        # A was predicted to lose:
      } else {
        inc(cont_table[2,1])
      }
      # A lost:
    } else {
      # A was predicted to win:
      if (identical(as.character(df$Species.A[i]), as.character(df[i, paste(law, ".Law.Prediction", sep = "")])) == TRUE) {
        inc(cont_table[1,2])
        # A was predicted to lose:
      } else {
        inc(cont_table[2,2])
      }
    }
  }
  return(cont_table)
}

chisq.test(populate_cont_table(both_group, "Linear"))
chisq.test(populate_cont_table(both_group, "Square"))
chisq.test(populate_cont_table(both_group2, "Linear"))
chisq.test(populate_cont_table(both_group2, "Square"))

# To distinguish between the linear and square laws, we run another chi-squared test.
# This time, the contingency table has the following structure: [1,1] correctly
# predicted by both models, [1,2] correctly predicted by the linear law but not the
# square law, [2,1] correctly predicted by the square law but not the linear law, 
# [2,2] predicted incorrectly by both laws.
law_vs_law <- function(df) {
  cont_table <- matrix(0, nrow = 2, ncol = 2)
  for(i in 1:nrow(df)) {
    if (df$Linear.True[i] == 1 & df$Square.True[i] == 1) {
      inc(cont_table[1,1])
    } else if (df$Linear.True[i] == 1 & df$Square.True[i] == 0) {
      inc(cont_table[1,2])
    } else if (df$Linear.True[i] == 0 & df$Square.True[i] == 1) {
      inc(cont_table[2,1])
    } else {
      inc(cont_table[2,2])
    }
  }
  return(cont_table)
}

laws_mass <- law_vs_law(both_group)
laws_mass2_3 <- law_vs_law(both_group2)
chisq.test(laws_mass)
chisq.test(laws_mass2_3)

Response to Behav. Ecol. reviews

Fish size correlation

We discovered a discrepancy between the identification of the longnose butterflyfish in our fish size data (FishSize.csv) and Table 1 of our manuscript: the same species was identified as Forcipiger longirostris in the former and as the closely related F. flavissimus in the latter. Based on the photograph we took of the species, we concluded that the latter identification was correct, as the lower half of its eye was silver rather than black (http://reefs.com/2015/12/31/long-nosed-butterflies-part-2-forcipiger/). We therefore replaced the maximum attained size reported for F. longirostris by Nelson (2005) with that reported for F. flavissimus. Below, we recalculate the correlation between the reported attained sizes and the observed sizes:

# Maximum attained sizes according to Randall (2005):
sizes <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", header = TRUE)
colnames(sizes)[1] <- "Species"
booksizes <- sizes$Attain[c(1:3,5:8,10:24)]

# Distributions of sizes observed in the field:
fishsize <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FieldFishSize.csv", sep = ",", strip.white = TRUE)

# Calculate the mean, median, and maximum of the observed size distributions:
meansizes <- colMeans(fishsize[,c(2:4,6:9,11:25)], na.rm = TRUE)
colMedians <- function(matrix) {
  return(apply(matrix, 2, median, na.rm = TRUE))
}
colMax <- function(matrix) {
  return(apply(matrix, 2, max, na.rm = TRUE))
}
mediansizes <- colMedians(fishsize[,c(2:4,6:9,11:25)])
maxsizes <- colMax(fishsize[,c(2:4,6:9,11:25)])

# Do book sizes correlate with live observation sizes?
g <- cor.test(booksizes, mediansizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, mediansizes, method = "spearman", :
## Cannot compute exact p-value with ties
paste("correlation: ", signif(g$estimate[[1]],5), ", p-value: ", g$p.value, sep = "")
## [1] "correlation: 0.71489, p-value: 0.000184888348145873"
h <- cor.test(booksizes, meansizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, meansizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with ties
paste("correlation: ", signif(h$estimate[[1]],5), ", p-value: ", h$p.value, sep = "")
## [1] "correlation: 0.76279, p-value: 3.66099985642471e-05"
i <- cor.test(booksizes, maxsizes, method = "spearman", alternative = "two.sided")
## Warning in cor.test.default(booksizes, maxsizes, method = "spearman",
## alternative = "two.sided"): Cannot compute exact p-value with ties
paste("correlation: ", signif(i$estimate[[1]],5), ", p-value: ", i$p.value, sep = "")
## [1] "correlation: 0.81934, p-value: 3.09742996782284e-06"
j <- cor.test(booksizes, mediansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(j$estimate[[1]],5), ", p-value: ", j$p.value, sep = "")
## [1] "correlation: 0.83931, p-value: 1.04831251950256e-06"
k <- cor.test(booksizes, meansizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(k$estimate[[1]],5), ", p-value: ", k$p.value, sep = "")
## [1] "correlation: 0.85355, p-value: 4.41082818979503e-07"
l <- cor.test(booksizes, maxsizes, method = "pearson", alternative = "two.sided")
paste("correlation: ", signif(l$estimate[[1]],5), ", p-value: ", l$p.value, sep = "")
## [1] "correlation: 0.83313, p-value: 1.48832376900543e-06"

CBI Calculation

Going back to Clutton-Brock et al. (1979), we discovered that our original implementation of the CBI-calculating function had an error. When counting the species beaten by those species that were dominated by the subject, we need to exclude the subject itself; the same is true when tallying the species that dominated those species which dominated the subject. This was not done in the original code, but it is done below:

# Remember:
# The A column represents the species that A dominates (A wins)
# The A row represents the species that dominate A (A losses)
# (B + b + 1)/(L + l + 1)
cbi_new <- function(df, species) {
  # number of entries in column 'species' that are neither 0 nor NA:
  B <- sum(df[,species] != 0, na.rm = TRUE)
  # species corresponding to those entries:
  dominates <- rownames(df)[which(df[,species] != 0)]
  # number of entries in the columns corresponding to those species that are neither 0 nor NA:
  if (length(dominates) > 0) {
    # exclude 'species' itself from consideration
    b <- sum(df[rownames(df) != species, colnames(df) %in% species] != 0, na.rm = TRUE)
  } else {
    b <- 0
  }
  # number of entries in row B that are neither 0 nor NA:
  L <- sum(df[species,], na.rm = TRUE)
  # species corresponding to those entries:
  dominated <- colnames(df)[which(df[species,] != 0)]
  # number of entries in the rows corresponding to those species that are neither 0 nor NA:
  if (length(dominated) > 0) {
    # exclude 'species' itself from consideration
    l <- sum(df[rownames(df) %in% dominated, colnames(df) != species] != 0, na.rm = TRUE)
  } else {
    l <- 0
  }
  return((B + b + 1)/(L + l + 1))
}

# Apply it to the full matrix:
cbi_hierarchy_new <- function(df) {
  cbi_matrix <- data.frame(matrix(NA, nrow = nrow(df), ncol = 2))
  colnames(cbi_matrix) <- c("Species", "CBI")
  for (i in 1:nrow(df)) {
    cbi_matrix$Species[i] <- rownames(df)[i]
    cbi_matrix$CBI[i] <- cbi_new(df, rownames(df)[i])
  }
  cbi_matrix <- cbi_matrix[order(cbi_matrix$CBI, decreasing = TRUE),]
  return(cbi_matrix)
}

Phylogenetic GLS

First, we need to drop the “Juvenile Parrotfish” category from the analysis, since it is a mixture of multiple species:

nojp_domall <- exclude(nocornet_domall, "Juvenile Parrotfish")
nojp_scoredall <- exclude(scoredall, "Juvenile Parrotfish")
nojp_bothall <- exclude(nocornet_bothall, "Juvenile Parrotfish")
nojp_live <- exclude(nocornet_live, "Juvenile Parrotfish")
nojp_scored <- exclude(scoredoooall, "Juvenile Parrotfish")
nojp_both <- exclude(nocornet_both, "Juvenile Parrotfish")

# Recalculate CBI:
cbi_nojp_domall <- cbi_hierarchy_new(nojp_domall)
cbi_nojp_scoredall <- cbi_hierarchy_new(nojp_scoredall)
cbi_nojp_bothall <- cbi_hierarchy_new(nojp_bothall)
cbi_nojp_ooolive <- cbi_hierarchy_new(nojp_live)
cbi_nojp_ooovideo <- cbi_hierarchy_new(nojp_scored)
cbi_nojp_oooboth <- cbi_hierarchy_new(nojp_both)

# Recalculate the frequency of winning and combine the results with the CBI hierarchy
# and mass-length data into a single data frame for further testing. We will again focus
# on one-on-one data only from now on.

# This is a generalized version of the no_surgeon() function defined above:
drop_taxa <- function(df, taxa_to_drop) {
  return(df[!df$Species.A %in% taxa_to_drop & !df$Species.B %in% taxa_to_drop,])
}

sheet_ooo5 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
scored_ooo5 <- drop_taxa(scored_ooo4, c("Cornetfish", "Juvenile Parrotfish"))
both_ooo5 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish"))

sheet_phylo_test <- merge(cbi_nojp_ooolive, win_freq(sheet_ooo5), by = "Species")
sheet_phylo_test <- merge(sheet_phylo_test, lengthmass, by = "Species")
# We need to provide row names to ensure that each row of the data frame is linked to the
# correct tree tip:
rownames(sheet_phylo_test) <- sheet_phylo_test$Species
scored_phylo_test <- merge(cbi_nojp_ooovideo, win_freq(scored_ooo5), by = "Species")
scored_phylo_test <- merge(scored_phylo_test, lengthmass, by = "Species")
rownames(scored_phylo_test) <- scored_phylo_test$Species
both_phylo_test <- merge(cbi_nojp_oooboth, win_freq(both_ooo5), by = "Species")
both_phylo_test <- merge(both_phylo_test, lengthmass, by = "Species")
rownames(both_phylo_test) <- both_phylo_test$Species

# Finally, repeat all of the above without the spotted porcupinefish (we will want to see
# how its exclusion as an outlier is going to change the results). This species was never
# observed on video, so there is no need to remove it from the video-only data. Note that
# the exclude() function is not safeguarded against attempts to remove taxa that are not
# present in the data frame to begin with.

nojpsp_domall <- exclude(nojp_domall, "Spotted Porcupinefish")
nojpsp_bothall <- exclude(nojp_bothall, "Spotted Porcupinefish")
nojpsp_live <- exclude(nojp_live, "Spotted Porcupinefish")
nojpsp_both <- exclude(nojp_both, "Spotted Porcupinefish")

cbi_nojpsp_domall <- cbi_hierarchy_new(nojpsp_domall)
cbi_nojpsp_bothall <- cbi_hierarchy_new(nojpsp_bothall)
cbi_nojpsp_ooolive <- cbi_hierarchy_new(nojpsp_live)
cbi_nojpsp_oooboth <- cbi_hierarchy_new(nojpsp_both)

sheet_ooo6 <- drop_taxa(sheet_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))
both_ooo6 <- drop_taxa(both_ooo4, c("Cornetfish", "Juvenile Parrotfish", "Spotted Porcupinefish"))

sheet_phylo_test_nosp <- merge(cbi_nojpsp_ooolive, win_freq(sheet_ooo6), by = "Species")
sheet_phylo_test_nosp <- merge(sheet_phylo_test_nosp, lengthmass, by = "Species")
rownames(sheet_phylo_test_nosp) <- sheet_phylo_test_nosp$Species
both_phylo_test_nosp <- merge(cbi_nojpsp_oooboth, win_freq(both_ooo6), by = "Species")
both_phylo_test_nosp <- merge(both_phylo_test_nosp, lengthmass, by = "Species")
rownames(both_phylo_test_nosp) <- both_phylo_test_nosp$Species

Recalculate the properties of the one-on-one matrices created by the exclusion of juvenile parrotfish:

paste("The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_live), "%", sep = "")
## [1] "The completeness of the live-only one-on-one matrix after the exclusion of juvenile parrotfish is 35.7142857142857%"
paste("The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_scored), "%", sep = "")
## [1] "The completeness of the video-only one-on-one matrix after the exclusion of juvenile parrotfish is 26.6666666666667%"
paste("The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is ", compl(nojp_both), "%", sep = "")
## [1] "The completeness of the combined one-on-one matrix after the exclusion of juvenile parrotfish is 34.6320346320346%"
devries(nojp_live, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3125917 
##  p-value= 0.0042

devries(nojp_scored, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.3038453 
##  p-value= 0.07

devries(nojp_both, Nperms = 10000, history = FALSE, plot = TRUE)
## h-modified = 0.2824597 
##  p-value= 0.0062

Now we can run phylogenetic regression:

library(nlme)
library(phytools)
## Loading required package: ape
## Loading required package: maps
fishtree <- read.tree("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/small_tree.tre")

# We will have to match the scientific names in the tip labels against the common names
# we use for the remaining data:
species_table <- read.csv("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/FishSize.csv", sep = ";", stringsAsFactors = FALSE)
fishtree$tip.label <- species_table$Common.Name[match(unlist(lapply(strsplit(fishtree$tip.label, "_"), function(x) paste(x[1], x[2], sep = " "))), species_table$Scientific.Name)]

# Make a new tree that does not include the spotted porcupinefish:
fishtree_nosp <- drop.tip(fishtree, "Spotted Porcupinefish")

# Run phylogenetic GLS on all parameter combinations (mass versus mass raised to two thirds, CBI vs. frequency of winning, spotted porcupinefish vs. no spotted porcupinefish):
pgls_CBI_Mass <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass)
## Generalized least squares fit by maximum likelihood
##   Model: CBI ~ Mass 
##   Data: both_phylo_test 
##        AIC      BIC   logLik
##   42.96499 47.32916 -17.4825
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##   lambda 
## 1.026069 
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept)  0.6181935 0.3625948  1.7049155  0.1037
## Mass        -0.0002343 0.0006324 -0.3705309  0.7149
## 
##  Correlation: 
##      (Intr)
## Mass -0.341
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.7446376 -0.6191146 -0.3996205 -0.2396558  3.0142764 
## 
## Residual standard error: 0.8129451 
## Degrees of freedom: 22 total; 20 residual
pgls_CBI_Mass2_3 <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_CBI_Mass2_3)
## Generalized least squares fit by maximum likelihood
##   Model: CBI ~ Mass2_3 
##   Data: both_phylo_test 
##        AIC     BIC    logLik
##   43.09373 47.4579 -17.54687
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##   lambda 
## 1.024889 
## 
## Coefficients:
##                  Value Std.Error    t-value p-value
## (Intercept)  0.5973099 0.3862164  1.5465678  0.1376
## Mass2_3     -0.0008940 0.0065095 -0.1373374  0.8921
## 
##  Correlation: 
##         (Intr)
## Mass2_3 -0.47 
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.7163492 -0.5910616 -0.4011444 -0.2178379  3.0026720 
## 
## Residual standard error: 0.8133159 
## Degrees of freedom: 22 total; 20 residual
pgls_WinFreq_Mass <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass)
## Generalized least squares fit by maximum likelihood
##   Model: WinFreq ~ Mass 
##   Data: both_phylo_test 
##        AIC      BIC     logLik
##   9.610968 13.97514 -0.8054841
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7608329 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept) 0.3867176 0.12406998 3.1169316  0.0054
## Mass        0.0001547 0.00024972 0.6193854  0.5427
## 
##  Correlation: 
##      (Intr)
## Mass -0.389
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3855278 -0.9836204 -0.3401884  0.1149432  1.9053666 
## 
## Residual standard error: 0.2980929 
## Degrees of freedom: 22 total; 20 residual
pgls_WinFreq_Mass2_3 <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree), data = both_phylo_test, method = "ML")
summary(pgls_WinFreq_Mass2_3)
## Generalized least squares fit by maximum likelihood
##   Model: WinFreq ~ Mass2_3 
##   Data: both_phylo_test 
##        AIC      BIC     logLik
##   9.341281 13.70545 -0.6706403
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7418624 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept) 0.3558832 0.13389644 2.6578989  0.0151
## Mass2_3     0.0021802 0.00266542 0.8179602  0.4230
## 
##  Correlation: 
##         (Intr)
## Mass2_3 -0.551
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4406599 -0.9320170 -0.4012155  0.1434963  1.9720109 
## 
## Residual standard error: 0.2934742 
## Degrees of freedom: 22 total; 20 residual
pgls_CBI_Mass_nosp <- gls(CBI ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass_nosp)
## Generalized least squares fit by maximum likelihood
##   Model: CBI ~ Mass 
##   Data: both_phylo_test_nosp 
##        AIC      BIC    logLik
##   38.77226 42.95035 -15.38613
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.9780776 
## 
## Coefficients:
##                 Value Std.Error  t-value p-value
## (Intercept) 0.4131129 0.3416725 1.209090  0.2415
## Mass        0.0017121 0.0013866 1.234727  0.2320
## 
##  Correlation: 
##      (Intr)
## Mass -0.476
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.9737380 -0.6492455 -0.4673056 -0.2472223  2.9127563 
## 
## Residual standard error: 0.7159469 
## Degrees of freedom: 21 total; 19 residual
pgls_CBI_Mass2_3_nosp <- gls(CBI ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_CBI_Mass2_3_nosp)
## Generalized least squares fit by maximum likelihood
##   Model: CBI ~ Mass2_3 
##   Data: both_phylo_test_nosp 
##        AIC      BIC    logLik
##   39.20199 43.38008 -15.60099
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.9938416 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 0.4139249 0.3738005 1.1073418  0.2820
## Mass2_3     0.0093035 0.0095417 0.9750351  0.3418
## 
##  Correlation: 
##         (Intr)
## Mass2_3 -0.551
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -0.8003119 -0.6502855 -0.5034040 -0.2894925  2.9422145 
## 
## Residual standard error: 0.7394883 
## Degrees of freedom: 21 total; 19 residual
pgls_WinFreq_Mass_nosp <- gls(WinFreq ~ Mass, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass_nosp)
## Generalized least squares fit by maximum likelihood
##   Model: WinFreq ~ Mass 
##   Data: both_phylo_test_nosp 
##        AIC      BIC     logLik
##   9.614419 13.79251 -0.8072095
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7513216 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept) 0.3365943 0.13968518 2.4096639  0.0263
## Mass        0.0006368 0.00066595 0.9561862  0.3510
## 
##  Correlation: 
##      (Intr)
## Mass -0.554
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4891060 -0.9730536 -0.5088269  0.0878157  2.0020512 
## 
## Residual standard error: 0.2987457 
## Degrees of freedom: 21 total; 19 residual
pgls_WinFreq_Mass2_3_nosp <- gls(WinFreq ~ Mass2_3, correlation = corPagel(1, phy = fishtree_nosp), data = both_phylo_test_nosp, method = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)
## Generalized least squares fit by maximum likelihood
##   Model: WinFreq ~ Mass2_3 
##   Data: both_phylo_test_nosp 
##        AIC      BIC     logLik
##   9.600858 13.77895 -0.8004288
## 
## Correlation Structure: corPagel
##  Formula: ~1 
##  Parameter estimate(s):
##    lambda 
## 0.7461108 
## 
## Coefficients:
##                  Value  Std.Error   t-value p-value
## (Intercept) 0.31370935 0.15285260 2.0523652  0.0542
## Mass2_3     0.00450538 0.00465324 0.9682252  0.3451
## 
##  Correlation: 
##         (Intr)
## Mass2_3 -0.653
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.5174832 -1.0098362 -0.5279994  0.0544028  1.9944100 
## 
## Residual standard error: 0.2978518 
## Degrees of freedom: 21 total; 19 residual

We can see that for all models involving CBI, we get very weak relationships (characterized by nearly horizontal regression lines) that are nevertheless extremely significant. A more careful examination reveals that this is because Pagel’s \(\lambda\) estimated for these models exceeds 1 (moreover, it is identical in all the four cases), even though it is only defined between 0 and 1. The cause of this is unclear (but likely related to the shape of the tree, especially to the presence of terminal branches longer than expected under the Yule process), as is the interpretability of the resulting estimate. Based on Liam Revell’s remarks in the first of the two links given below, we turn to the caper package to fit a phylogenetic GLS under the constraint that \(\lambda \leq 1\):

http://blog.phytools.org/2016/01/phylogenetic-regression-when-estimated.html http://grokbase.com/t/r/r-sig-phylo/113s5mhjgg/pagels-lambda-greater-than-one

if(!require("caper")) {install.packages("caper")}
## Loading required package: caper
## Loading required package: MASS
## Loading required package: mvtnorm
library(caper)
both_comparative <- comparative.data(data = both_phylo_test, phy = fishtree, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
both_comparative_nosp <- comparative.data(data = both_phylo_test_nosp, phy = fishtree_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)

pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass)
## 
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14667 -0.02814  0.01556  0.05592  0.08199 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.010507
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.524, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate  Std. Error t value Pr(>|t|)  
## (Intercept)  0.61306817  0.34503191  1.7768  0.09081 .
## Mass        -0.00020831  0.00061615 -0.3381  0.73882  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07567 on 20 degrees of freedom
## Multiple R-squared: 0.005683,    Adjusted R-squared: -0.04403 
## F-statistic: 0.1143 on 1 and 20 DF,  p-value: 0.7388
pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_CBI_Mass2_3)
## 
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14586 -0.02786  0.01345  0.05485  0.09016 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.016049
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.453, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate  Std. Error t value Pr(>|t|)
## (Intercept)  0.58731317  0.37068156  1.5844   0.1288
## Mass2_3     -0.00053638  0.00644425 -0.0832   0.9345
## 
## Residual standard error: 0.07587 on 20 degrees of freedom
## Multiple R-squared: 0.0003463,   Adjusted R-squared: -0.04964 
## F-statistic: 0.006928 on 1 and 20 DF,  p-value: 0.9345
pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass_nosp)
## 
## Call:
## pgls(formula = CBI ~ Mass, data = both_comparative_nosp, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14830 -0.04944 -0.01476  0.01939  0.18616 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.978
##    lower bound : 0.000, p = 0.0081072
##    upper bound : 1.000, p = 0.83889
##    95.0% CI   : (0.458, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4131104  0.3416712  1.2091   0.2415
## Mass        0.0017121  0.0013866  1.2347   0.2320
## 
## Residual standard error: 0.06986 on 19 degrees of freedom
## Multiple R-squared: 0.07428, Adjusted R-squared: 0.02556 
## F-statistic: 1.525 on 1 and 19 DF,  p-value: 0.232
pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_CBI_Mass2_3_nosp)
## 
## Call:
## pgls(formula = CBI ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.151213 -0.018977  0.009497  0.064289  0.139694 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.994
##    lower bound : 0.000, p = 0.0090652
##    upper bound : 1.000, p = 0.95565
##    95.0% CI   : (0.461, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4139220  0.3737994  1.1073   0.2820
## Mass2_3     0.0093036  0.0095417  0.9750   0.3418
## 
## Residual standard error: 0.07216 on 19 degrees of freedom
## Multiple R-squared: 0.04765, Adjusted R-squared: -0.00247 
## F-statistic: 0.9507 on 1 and 19 DF,  p-value: 0.3418
pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass)
## 
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.067829 -0.016253  0.002191  0.017853  0.053835 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.761
##    lower bound : 0.000, p = 0.026067
##    upper bound : 1.000, p = 0.060639
##    95.0% CI   : (0.083, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 0.38671755 0.12406983  3.1169 0.005432 **
## Mass        0.00015467 0.00024972  0.6194 0.542654   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02902 on 20 degrees of freedom
## Multiple R-squared: 0.01882, Adjusted R-squared: -0.03024 
## F-statistic: 0.3836 on 1 and 20 DF,  p-value: 0.5427
pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
summary(pgls_WinFreq_Mass2_3)
## 
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.056556 -0.021804 -0.006782  0.014359  0.038175 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.742
##    lower bound : 0.000, p = 0.028143
##    upper bound : 1.000, p = 0.052267
##    95.0% CI   : (0.071, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.3558831  0.1338963  2.6579   0.0151 *
## Mass2_3     0.0021802  0.0026654  0.8180   0.4230  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02857 on 20 degrees of freedom
## Multiple R-squared: 0.03237, Adjusted R-squared: -0.01601 
## F-statistic: 0.6691 on 1 and 20 DF,  p-value: 0.423
pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass_nosp)
## 
## Call:
## pgls(formula = WinFreq ~ Mass, data = both_comparative_nosp, 
##     lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.063430 -0.018323 -0.003256  0.012007  0.054953 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.751
##    lower bound : 0.000, p = 0.022585
##    upper bound : 1.000, p = 0.055727
##    95.0% CI   : (0.103, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.33659421 0.13968505  2.4097  0.02627 *
## Mass        0.00063677 0.00066595  0.9562  0.35099  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02915 on 19 degrees of freedom
## Multiple R-squared: 0.04591, Adjusted R-squared: -0.004304 
## F-statistic: 0.9143 on 1 and 19 DF,  p-value: 0.351
pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp, lambda = "ML")
summary(pgls_WinFreq_Mass2_3_nosp)
## 
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = both_comparative_nosp, 
##     lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.061849 -0.008427  0.005194  0.019045  0.054724 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.746
##    lower bound : 0.000, p = 0.023835
##    upper bound : 1.000, p = 0.055536
##    95.0% CI   : (0.095, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.3137092  0.1528525  2.0524  0.05417 .
## Mass2_3     0.0045054  0.0046532  0.9682  0.34510  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02906 on 19 degrees of freedom
## Multiple R-squared: 0.04702, Adjusted R-squared: -0.003137 
## F-statistic: 0.9375 on 1 and 19 DF,  p-value: 0.3451
# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(pgls_CBI_Mass)

par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3)

par(mfrow=c(2,2))
plot(pgls_CBI_Mass_nosp)

par(mfrow=c(2,2))
plot(pgls_CBI_Mass2_3_nosp)

par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass)

par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3)

par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass_nosp)

par(mfrow=c(2,2))
plot(pgls_WinFreq_Mass2_3_nosp)

This worked great: the \(\lambda\) and p values for the winning frequency are the same as before, and while \(\lambda\) is still running up against the upper bound of 1 for CBI, the fact that it cannot exceed it results in much more sensible p values. We can now proceed to visualizing the results:

# Generate a new Figure 1:

png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Figure_1.png", width = 3600, height = 3600, pointsize = 96)
# setEPS()
# postscript("/Users/David/Documents/College/2017–18/Winter 2018/FBQ_git/Figure_1.eps")
spec_dec <- function(x, k) trimws(format(round(x, k), nsmall = k)) # credit: Jeromy Anglim, Stack Overflow
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))

par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(pgls_CBI_Mass, lwd = 5)
abline(pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass$param[2]
pseudor2 <- summary(pgls_CBI_Mass)$r.squared
pvalue <- summary(pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.3, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.1, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.85, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_CBI_Mass2_3, lwd = 5)
abline(pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.8, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 1.55, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass, both_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass, lwd = 5)
abline(pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass)$r.squared
pvalue <- summary(pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.7, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.7 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.7 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(both_phylo_test$Mass2_3, both_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(both_phylo_test$Mass == both_phylo_test$Mass[both_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(pgls_WinFreq_Mass2_3, lwd = 5)
abline(pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.98, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.98 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.98 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.225, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.225 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()

Group size-only and body size-only models

Below, we conduct another test in which the \(\alpha\) coefficients (individual fighting abilities) are all set to 1, and therefore effectively disregarded. We also evaluate the fit of another model, which includes body size as the only parameter. Note that neither model requires us to repeat the test for body mass raised to 2/3. The group size-only model does not include this parameter at all, while the body size-only model will always yield the same predictions under both mass scalings, since \(m^{\frac{2}{3}}\) is a monotonically increasing function of \(m\).

# Add a new column to both_group with predictions based on group size alone:
for(i in 1:nrow(both_group)) {
  if (both_group$Group.A[i] > both_group$Group.B[i]) {
    both_group[i,14] <- both_group$Species.A[i] 
  } else {
    both_group[i,14] <- both_group$Species.B[i] 
  }
}

# Add a new column to both_group with predictions based on body size alone:
for(i in 1:nrow(both_group)) {
  if (both_group$BodySize.A[i] > both_group$BodySize.B[i]) {
    both_group[i,15] <- both_group$Species.A[i] 
  } else {
    both_group[i,15] <- both_group$Species.B[i] 
  }
}
colnames(both_group)[14:15] <- c("Group.Size.Only.Prediction","Body.Size.Only.Prediction")

# Test if the group size-only model works:
for(i in 1:nrow(both_group)) {
  if (identical(as.character(both_group$Winner[i]), as.character(both_group$Group.Size.Only.Prediction[i])) == TRUE) {
    both_group[i,16] <- 1
  } else {
    both_group[i,16] <- 0
  }
}

# Test if the body size-only model works:
for(i in 1:nrow(both_group)) {
  if (identical(as.character(both_group$Winner[i]), as.character(both_group$Body.Size.Only.Prediction[i])) == TRUE) {
    both_group[i,17] <- 1
  } else {
    both_group[i,17] <- 0
  }
}
colnames(both_group)[16:17] <- c("Group.Size.Only.True", "Body.Size.Only.True")

# Proportion of successful predictions:
sum(both_group$Group.Size.Only.True)/nrow(both_group)
## [1] 0.2289157
sum(both_group$Body.Size.Only.True)/nrow(both_group)
## [1] 0.6987952
# Test significance:
dbinom(sum(both_group$Group.Size.Only.True), size = nrow(both_group), prob = 1/2)
## [1] 2.643039e-07
dbinom(sum(both_group$Body.Size.Only.True), size = nrow(both_group), prob = 1/2)
## [1] 0.0001118917

In what cases did the body size-only model work but the next best model (Lanchester’s linear law) did not?

both_group[both_group$Body.Size.Only.True == 1 & both_group$Linear.True == 0,]
##                    Species.A             Species.B Group.A Group.B
## 59                Damselfish  Scissortail Sergeant       1       2
## 87                Damselfish         Sixbar Wrasse       1       2
## 246               Damselfish   Juvenile Parrotfish       1      20
## 259          Regal Angelfish            Brown Tang       1       2
## 273 Orange-lined Triggerfish Threeband Pennantfish       1       2
## 276 Orange-lined Triggerfish            Brown Tang       1       2
## 324               Damselfish  Scissortail Sergeant       1       3
## 329               Damselfish         Sixbar Wrasse       1       2
## 487               Damselfish   Juvenile Parrotfish       1       7
## 515               Damselfish   Juvenile Parrotfish       1      14
## 554               Damselfish  Scissortail Sergeant       1       2
##                       Winner BodySize.A BodySize.B Linear.Law.Prediction
## 59                Damselfish   38.88513   25.35153  Scissortail Sergeant
## 87                Damselfish   38.88513   22.72689         Sixbar Wrasse
## 246               Damselfish   38.88513    7.49326   Juvenile Parrotfish
## 259          Regal Angelfish  128.17052   78.11946            Brown Tang
## 273 Orange-lined Triggerfish  108.96596   87.65761 Threeband Pennantfish
## 276 Orange-lined Triggerfish  108.96596   78.11946            Brown Tang
## 324               Damselfish   38.88513   25.35153  Scissortail Sergeant
## 329               Damselfish   38.88513   22.72689         Sixbar Wrasse
## 487               Damselfish   38.88513    7.49326   Juvenile Parrotfish
## 515               Damselfish   38.88513    7.49326   Juvenile Parrotfish
## 554               Damselfish   38.88513   25.35153  Scissortail Sergeant
##     Square.Law.Prediction Linear.True Square.True
## 59   Scissortail Sergeant           0           0
## 87          Sixbar Wrasse           0           0
## 246   Juvenile Parrotfish           0           0
## 259            Brown Tang           0           0
## 273 Threeband Pennantfish           0           0
## 276            Brown Tang           0           0
## 324  Scissortail Sergeant           0           0
## 329         Sixbar Wrasse           0           0
## 487   Juvenile Parrotfish           0           0
## 515   Juvenile Parrotfish           0           0
## 554  Scissortail Sergeant           0           0
##                            Diet.A                        Diet.B
## 59                       Omnivore                      Omnivore
## 87                       Omnivore Benthic Invertebrate Consumer
## 246                      Omnivore         Herbivore/Detritivore
## 259 Benthic Invertebrate Consumer         Herbivore/Detritivore
## 273                      Omnivore                   Corallivore
## 276                      Omnivore         Herbivore/Detritivore
## 324                      Omnivore                      Omnivore
## 329                      Omnivore Benthic Invertebrate Consumer
## 487                      Omnivore         Herbivore/Detritivore
## 515                      Omnivore         Herbivore/Detritivore
## 554                      Omnivore                      Omnivore
##     Group.Size.Only.Prediction Body.Size.Only.Prediction
## 59        Scissortail Sergeant                Damselfish
## 87               Sixbar Wrasse                Damselfish
## 246        Juvenile Parrotfish                Damselfish
## 259                 Brown Tang           Regal Angelfish
## 273      Threeband Pennantfish  Orange-lined Triggerfish
## 276                 Brown Tang  Orange-lined Triggerfish
## 324       Scissortail Sergeant                Damselfish
## 329              Sixbar Wrasse                Damselfish
## 487        Juvenile Parrotfish                Damselfish
## 515        Juvenile Parrotfish                Damselfish
## 554       Scissortail Sergeant                Damselfish
##     Group.Size.Only.True Body.Size.Only.True
## 59                     0                   1
## 87                     0                   1
## 246                    0                   1
## 259                    0                   1
## 273                    0                   1
## 276                    0                   1
## 324                    0                   1
## 329                    0                   1
## 487                    0                   1
## 515                    0                   1
## 554                    0                   1

8 out of the 11 interactions in this category involved damselfish. What happens if we remove them from the dataset?

both_group3 <- both_group[both_group$Species.A != "Damselfish" & both_group$Species.B != "Damselfish", ]
nrow(both_group3)
## [1] 43
# Proportion of successful predictions:
sum(both_group3$Linear.True)/nrow(both_group3)
## [1] 0.8372093
sum(both_group3$Square.True)/nrow(both_group3)
## [1] 0.7209302
sum(both_group3$Group.Size.Only.True)/nrow(both_group3)
## [1] 0.255814
sum(both_group3$Body.Size.Only.True)/nrow(both_group3)
## [1] 0.7906977
# Test significance:
dbinom(sum(both_group3$Linear.True), size = nrow(both_group3), prob = 1/2)
## [1] 3.663458e-06
dbinom(sum(both_group3$Square.True), size = nrow(both_group3), prob = 1/2)
## [1] 0.001743806
dbinom(sum(both_group3$Group.Size.Only.True), size = nrow(both_group3), prob = 1/2)
## [1] 0.0006539272
dbinom(sum(both_group3$Body.Size.Only.True), size = nrow(both_group3), prob = 1/2)
## [1] 6.411051e-05
# With mass raised to the two-thirds power. Remember, this can change the results for the
# linear and square laws, but not for the body size-only or group size-only models.
both_group4 <- both_group2[both_group2$Species.A != "Damselfish" & both_group2$Species.B != "Damselfish" ,]
nrow(both_group4)
## [1] 43
# Proportion of successful predictions:
sum(both_group4$Linear.True)/nrow(both_group4)
## [1] 0.7906977
sum(both_group4$Square.True)/nrow(both_group4)
## [1] 0.5348837
# Test significance:
dbinom(sum(both_group4$Linear.True), size = nrow(both_group4), prob = 1/2)
## [1] 6.411051e-05
dbinom(sum(both_group4$Square.True), size = nrow(both_group4), prob = 1/2)
## [1] 0.1092038

In the next step, we calculate the statistical power of the linear law, square law, and body size-only model tests. All of the below is based on http://www.calvin.edu/~scofield/courses/m343/F15/handouts/binomialTestPower.pdf

# First, find the rejection region for n = 83 (the number of data points in both_group)
# under the null assumption (probability of successful prediction = 0.5) and significance
# level alpha = 0.05:
qbinom(0.025, nrow(both_group), 0.5) # 33
## [1] 33
qbinom(0.975, nrow(both_group), 0.5) # 50
## [1] 50
# Since P(33 \leq x \geq 50) is probably not equal to exactly 0.95, we need to find out 
# if the sum of the values below is larger or smaller than this significance threshold:
pbinom(33, nrow(both_group), 0.5) + (1 - pbinom(50, nrow(both_group), 0.5))
## [1] 0.06297226
pbinom(32, nrow(both_group), 0.5) + (1 - pbinom(51, nrow(both_group), 0.5))
## [1] 0.03752953
# Let's use 32 and 51 as the boundaries of the rejection region.

# Compute power as 1 - beta:
1 - (pbinom(51, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)) -
       pbinom(32, nrow(both_group), sum(both_group$Body.Size.Only.True)/nrow(both_group)))
## [1] 0.9377759

Compute the maximum likelihood estimate under the binomial and use it to compute AIC scores:

# -2*log(L^hat) for the body size-only model:
-2*log(dbinom(sum(both_group$Body.Size.Only.True), 
              nrow(both_group),
              sum(both_group$Body.Size.Only.True)/nrow(both_group)))
## [1] 4.705887
# -2*log(L^hat) for the linear law:
-2*log(dbinom(sum(both_group$Linear.True), 
              nrow(both_group),
              sum(both_group$Linear.True)/nrow(both_group)))
## [1] 4.765549

Finally, see how the body size-only and group-size only models perform when fitted to the datasets subsampled by diet:

test_omni
##                    Species.A            Species.B Group.A Group.B
## 59                Damselfish Scissortail Sergeant       1       2
## 113             Moorish Idol  Lemonpeel Angelfish       1       2
## 132               Damselfish  Lemonpeel Angelfish       1       2
## 187          Dusky Angelfish Scissortail Sergeant       1       3
## 288 Orange-lined Triggerfish  Lemonpeel Angelfish       1       2
## 324               Damselfish Scissortail Sergeant       1       3
## 353 Orange-lined Triggerfish  Lemonpeel Angelfish       1       3
## 355               Damselfish  Lemonpeel Angelfish       1       3
## 357               Damselfish  Lemonpeel Angelfish       1       3
## 358               Damselfish  Lemonpeel Angelfish       2       1
## 363 Orange-lined Triggerfish  Lemonpeel Angelfish       1       3
## 366             Moorish Idol  Lemonpeel Angelfish       1       3
## 368               Damselfish  Lemonpeel Angelfish       1       3
## 373               Damselfish  Lemonpeel Angelfish       1       3
## 440               Damselfish  Lemonpeel Angelfish       1       2
## 554               Damselfish Scissortail Sergeant       1       2
##                       Winner BodySize.A BodySize.B
## 59                Damselfish   38.88513  25.351526
## 113             Moorish Idol  171.38610   7.490788
## 132               Damselfish   38.88513   7.490788
## 187     Scissortail Sergeant   10.99104  25.351526
## 288 Orange-lined Triggerfish  108.96596   7.490788
## 324               Damselfish   38.88513  25.351526
## 353 Orange-lined Triggerfish  108.96596   7.490788
## 355               Damselfish   38.88513   7.490788
## 357               Damselfish   38.88513   7.490788
## 358               Damselfish   38.88513   7.490788
## 363 Orange-lined Triggerfish  108.96596   7.490788
## 366             Moorish Idol  171.38610   7.490788
## 368               Damselfish   38.88513   7.490788
## 373               Damselfish   38.88513   7.490788
## 440               Damselfish   38.88513   7.490788
## 554               Damselfish   38.88513  25.351526
##        Linear.Law.Prediction    Square.Law.Prediction Linear.True
## 59      Scissortail Sergeant     Scissortail Sergeant           0
## 113             Moorish Idol             Moorish Idol           1
## 132               Damselfish               Damselfish           1
## 187     Scissortail Sergeant     Scissortail Sergeant           1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish           1
## 324     Scissortail Sergeant     Scissortail Sergeant           0
## 353 Orange-lined Triggerfish Orange-lined Triggerfish           1
## 355               Damselfish      Lemonpeel Angelfish           1
## 357               Damselfish      Lemonpeel Angelfish           1
## 358               Damselfish               Damselfish           1
## 363 Orange-lined Triggerfish Orange-lined Triggerfish           1
## 366             Moorish Idol             Moorish Idol           1
## 368               Damselfish      Lemonpeel Angelfish           1
## 373               Damselfish      Lemonpeel Angelfish           1
## 440               Damselfish               Damselfish           1
## 554     Scissortail Sergeant     Scissortail Sergeant           0
##     Square.True   Diet.A   Diet.B
## 59            0 Omnivore Omnivore
## 113           1 Omnivore Omnivore
## 132           1 Omnivore Omnivore
## 187           1 Omnivore Omnivore
## 288           1 Omnivore Omnivore
## 324           0 Omnivore Omnivore
## 353           1 Omnivore Omnivore
## 355           0 Omnivore Omnivore
## 357           0 Omnivore Omnivore
## 358           1 Omnivore Omnivore
## 363           1 Omnivore Omnivore
## 366           1 Omnivore Omnivore
## 368           0 Omnivore Omnivore
## 373           0 Omnivore Omnivore
## 440           1 Omnivore Omnivore
## 554           0 Omnivore Omnivore
test_omni2
##                    Species.A            Species.B Group.A Group.B
## 59                Damselfish Scissortail Sergeant       1       2
## 113             Moorish Idol  Lemonpeel Angelfish       1       2
## 132               Damselfish  Lemonpeel Angelfish       1       2
## 187          Dusky Angelfish Scissortail Sergeant       1       3
## 288 Orange-lined Triggerfish  Lemonpeel Angelfish       1       2
## 324               Damselfish Scissortail Sergeant       1       3
## 353 Orange-lined Triggerfish  Lemonpeel Angelfish       1       3
## 355               Damselfish  Lemonpeel Angelfish       1       3
## 357               Damselfish  Lemonpeel Angelfish       1       3
## 358               Damselfish  Lemonpeel Angelfish       2       1
## 363 Orange-lined Triggerfish  Lemonpeel Angelfish       1       3
## 366             Moorish Idol  Lemonpeel Angelfish       1       3
## 368               Damselfish  Lemonpeel Angelfish       1       3
## 373               Damselfish  Lemonpeel Angelfish       1       3
## 440               Damselfish  Lemonpeel Angelfish       1       2
## 554               Damselfish Scissortail Sergeant       1       2
##                       Winner BodySize.A BodySize.B
## 59                Damselfish  11.477723   8.629840
## 113             Moorish Idol  30.854397   3.828409
## 132               Damselfish  11.477723   3.828409
## 187     Scissortail Sergeant   4.943401   8.629840
## 288 Orange-lined Triggerfish  22.813604   3.828409
## 324               Damselfish  11.477723   8.629840
## 353 Orange-lined Triggerfish  22.813604   3.828409
## 355               Damselfish  11.477723   3.828409
## 357               Damselfish  11.477723   3.828409
## 358               Damselfish  11.477723   3.828409
## 363 Orange-lined Triggerfish  22.813604   3.828409
## 366             Moorish Idol  30.854397   3.828409
## 368               Damselfish  11.477723   3.828409
## 373               Damselfish  11.477723   3.828409
## 440               Damselfish  11.477723   3.828409
## 554               Damselfish  11.477723   8.629840
##        Linear.Law.Prediction    Square.Law.Prediction Linear.True
## 59      Scissortail Sergeant     Scissortail Sergeant           0
## 113             Moorish Idol             Moorish Idol           1
## 132               Damselfish      Lemonpeel Angelfish           1
## 187     Scissortail Sergeant     Scissortail Sergeant           1
## 288 Orange-lined Triggerfish Orange-lined Triggerfish           1
## 324     Scissortail Sergeant     Scissortail Sergeant           0
## 353 Orange-lined Triggerfish      Lemonpeel Angelfish           1
## 355      Lemonpeel Angelfish      Lemonpeel Angelfish           0
## 357      Lemonpeel Angelfish      Lemonpeel Angelfish           0
## 358               Damselfish               Damselfish           1
## 363 Orange-lined Triggerfish      Lemonpeel Angelfish           1
## 366             Moorish Idol      Lemonpeel Angelfish           1
## 368      Lemonpeel Angelfish      Lemonpeel Angelfish           0
## 373      Lemonpeel Angelfish      Lemonpeel Angelfish           0
## 440               Damselfish      Lemonpeel Angelfish           1
## 554     Scissortail Sergeant     Scissortail Sergeant           0
##     Square.True   Diet.A   Diet.B
## 59            0 Omnivore Omnivore
## 113           1 Omnivore Omnivore
## 132           0 Omnivore Omnivore
## 187           1 Omnivore Omnivore
## 288           1 Omnivore Omnivore
## 324           0 Omnivore Omnivore
## 353           0 Omnivore Omnivore
## 355           0 Omnivore Omnivore
## 357           0 Omnivore Omnivore
## 358           1 Omnivore Omnivore
## 363           0 Omnivore Omnivore
## 366           0 Omnivore Omnivore
## 368           0 Omnivore Omnivore
## 373           0 Omnivore Omnivore
## 440           0 Omnivore Omnivore
## 554           0 Omnivore Omnivore
test_herb
##                  Species.A           Species.B Group.A Group.B
## 116             Brown Tang Juvenile Parrotfish       1       4
## 182 Dark-capped Parrotfish Juvenile Parrotfish       1      10
## 264  Bullethead Parrotfish          Brown Tang       1       2
## 282 Dark-capped Parrotfish Juvenile Parrotfish       1      20
##                     Winner BodySize.A BodySize.B  Linear.Law.Prediction
## 116             Brown Tang   78.11946    7.49326             Brown Tang
## 182 Dark-capped Parrotfish  311.91661    7.49326 Dark-capped Parrotfish
## 264  Bullethead Parrotfish  339.38694   78.11946  Bullethead Parrotfish
## 282 Dark-capped Parrotfish  311.91661    7.49326 Dark-capped Parrotfish
##     Square.Law.Prediction Linear.True Square.True                Diet.A
## 116   Juvenile Parrotfish           1           0 Herbivore/Detritivore
## 182   Juvenile Parrotfish           1           0 Herbivore/Detritivore
## 264 Bullethead Parrotfish           1           1 Herbivore/Detritivore
## 282   Juvenile Parrotfish           1           0 Herbivore/Detritivore
##                    Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivore
test_herb2
##                  Species.A           Species.B Group.A Group.B
## 116             Brown Tang Juvenile Parrotfish       1       4
## 182 Dark-capped Parrotfish Juvenile Parrotfish       1      10
## 264  Bullethead Parrotfish          Brown Tang       1       2
## 282 Dark-capped Parrotfish Juvenile Parrotfish       1      20
##                     Winner BodySize.A BodySize.B  Linear.Law.Prediction
## 116             Brown Tang   18.27425   3.829251             Brown Tang
## 182 Dark-capped Parrotfish   45.99306   3.829251 Dark-capped Parrotfish
## 264  Bullethead Parrotfish   48.65529  18.274246  Bullethead Parrotfish
## 282 Dark-capped Parrotfish   45.99306   3.829251    Juvenile Parrotfish
##     Square.Law.Prediction Linear.True Square.True                Diet.A
## 116   Juvenile Parrotfish           1           0 Herbivore/Detritivore
## 182   Juvenile Parrotfish           1           0 Herbivore/Detritivore
## 264            Brown Tang           1           0 Herbivore/Detritivore
## 282   Juvenile Parrotfish           0           0 Herbivore/Detritivore
##                    Diet.B
## 116 Herbivore/Detritivore
## 182 Herbivore/Detritivore
## 264 Herbivore/Detritivore
## 282 Herbivore/Detritivore
# The function below adds columns "Group.Size.Only.Prediction",
# "Body.Size.Only.Prediction", "Group.Size.Only.True", and "Body.Size.Only.True" to an 
# arbitrary data frame containing columns named "Group.A", "Group.B", "BodySizeA",
# "BodySizeB", and "Winner":
test_extra_two_models <- function(dataframe) {
  n <- ncol(dataframe)
  for(i in 1:nrow(dataframe)) {
    if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
      dataframe[i, n+1] <- dataframe$Species.A[i]
    } else {
      dataframe[i, n+1] <- dataframe$Species.B[i]
    }
    if (dataframe$BodySize.A[i] > dataframe$BodySize.B[i]) {
      dataframe[i, n+2] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+2] <- dataframe$Species.B[i] 
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+1])) == TRUE) {
      dataframe[i, n+3] <- 1
    } else {
      dataframe[i, n+3] <- 0
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+2])) == TRUE) {
      dataframe[i, n+4] <- 1
    } else {
      dataframe[i, n+4] <- 0
    }
  }
  colnames(dataframe)[(n+1):(n+4)] <- c("Group.Size.Only.Prediction",
                                        "Body.Size.Only.Prediction",
                                        "Group.Size.Only.True",
                                        "Body.Size.Only.True")
  return(dataframe)
}

test_omni <- test_extra_two_models(test_omni)
test_herb <- test_extra_two_models(test_herb)

# Proportion of successful predictions:
sum(test_omni$Group.Size.Only.True)/nrow(test_omni)
## [1] 0.125
sum(test_omni$Body.Size.Only.True)/nrow(test_omni)
## [1] 1
sum(test_herb$Group.Size.Only.True)/nrow(test_herb)
## [1] 0
sum(test_herb$Body.Size.Only.True)/nrow(test_herb)
## [1] 1
# Test significance:
dbinom(sum(test_omni$Group.Size.Only.True), size = nrow(test_omni), prob = 1/2)
## [1] 0.001831055
dbinom(sum(test_omni$Body.Size.Only.True), size = nrow(test_omni), prob = 1/2)
## [1] 1.525879e-05
dbinom(sum(test_herb$Group.Size.Only.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625
dbinom(sum(test_herb$Body.Size.Only.True), size = nrow(test_herb), prob = 1/2)
## [1] 0.0625

Live/video differences

We have already alculated the properties of matrices based on live-only and video-only observations. In the next step, we compare the hierarchies:

# Subsample each taxon-by-taxon dominance matrix down to the shared set of 15 species
nojp_live_shared <- exclude(nojp_live, setdiff(rownames(nojp_live), rownames(nojp_scored)))
nojp_scored_shared <- exclude(nojp_scored, setdiff(rownames(nojp_scored), rownames(nojp_live)))
nojp_both_shared <- exclude(nojp_both, rownames(nojp_both)[!rownames(nojp_both) %in% intersect(rownames(nojp_live), rownames(nojp_scored))])

# Recalculate the CBI for the 15 species based on the live-only and video-only datasets:
live_subsampled_cbi <- cbi_hierarchy_new(nojp_live_shared)
colnames(live_subsampled_cbi)[2] <- "Live"
video_subsampled_cbi <- cbi_hierarchy_new(nojp_scored_shared)
colnames(video_subsampled_cbi)[2] <- "Video"
subsampled_cbi <- merge(live_subsampled_cbi, video_subsampled_cbi, by = "Species")
both_subsampled_cbi <- cbi_hierarchy_new(nojp_both_shared)
subsampled_cbi <- merge(subsampled_cbi, both_subsampled_cbi, by = "Species")
colnames(subsampled_cbi)[4] <- "Both"

# Recalculate the frequency of winning for the 15 species based on the live-only and
# video-only datasets:
live_subsampled_winfreq <- win_freq(sheet_ooo5)[win_freq(sheet_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
live_subsampled_winfreq <- live_subsampled_winfreq[order(live_subsampled_winfreq$Species),]
colnames(live_subsampled_winfreq)[2] <- "Live"
video_subsampled_winfreq <- win_freq(scored_ooo5)[win_freq(scored_ooo5)$Species %in% intersect(win_freq(sheet_ooo5)$Species, win_freq(scored_ooo5)$Species),]
video_subsampled_winfreq <- video_subsampled_winfreq[order(video_subsampled_winfreq$Species),]
colnames(video_subsampled_winfreq)[2] <- "Video"
subsampled_winfreq <- merge(live_subsampled_winfreq, video_subsampled_winfreq, by = "Species")
subsampled_winfreq <- merge(subsampled_winfreq, win_freq(both_ooo5), by = "Species")
colnames(subsampled_winfreq)[4] <- "Both"

if(!require("RcmdrMisc")) {install.packages("RcmdrMisc")} # Beware: has many dependencies!
## Loading required package: RcmdrMisc
## Loading required package: car
## Loading required package: carData
## Loading required package: sandwich
library(RcmdrMisc)

rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "spearman")
## 
##  Spearman correlations:
##         Live  Video   Both
## Live  1.0000 0.2898 0.9357
## Video 0.2898 1.0000 0.3792
## Both  0.9357 0.3792 1.0000
## 
##  Number of observations: 15 
## 
##  Pairwise two-sided p-values:
##       Live   Video  Both  
## Live         0.2948 <.0001
## Video 0.2948        0.1633
## Both  <.0001 0.1633       
## 
##  Adjusted p-values (Holm's method)
##       Live   Video  Both  
## Live         0.3266 <.0001
## Video 0.3266        0.3266
## Both  <.0001 0.3266
rcorr.adjust(as.matrix(subsampled_cbi[,-1]), type = "pearson")
## 
##  Pearson correlations:
##         Live  Video   Both
## Live  1.0000 0.6703 0.8643
## Video 0.6703 1.0000 0.8914
## Both  0.8643 0.8914 1.0000
## 
##  Number of observations: 15 
## 
##  Pairwise two-sided p-values:
##       Live   Video  Both  
## Live         0.0062 <.0001
## Video 0.0062        <.0001
## Both  <.0001 <.0001       
## 
##  Adjusted p-values (Holm's method)
##       Live   Video  Both  
## Live         0.0062 <.0001
## Video 0.0062        <.0001
## Both  <.0001 <.0001
rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "spearman")
## 
##  Spearman correlations:
##         Live  Video   Both
## Live  1.0000 0.6188 0.8959
## Video 0.6188 1.0000 0.8307
## Both  0.8959 0.8307 1.0000
## 
##  Number of observations: 15 
## 
##  Pairwise two-sided p-values:
##       Live   Video  Both  
## Live         0.0139 <.0001
## Video 0.0139        0.0001
## Both  <.0001 0.0001       
## 
##  Adjusted p-values (Holm's method)
##       Live   Video  Both  
## Live         0.0139 <.0001
## Video 0.0139        0.0003
## Both  <.0001 0.0003
rcorr.adjust(as.matrix(subsampled_winfreq[,-1]), type = "pearson")
## 
##  Pearson correlations:
##         Live  Video   Both
## Live  1.0000 0.6208 0.9133
## Video 0.6208 1.0000 0.8731
## Both  0.9133 0.8731 1.0000
## 
##  Number of observations: 15 
## 
##  Pairwise two-sided p-values:
##       Live   Video  Both  
## Live         0.0135 <.0001
## Video 0.0135        <.0001
## Both  <.0001 <.0001       
## 
##  Adjusted p-values (Holm's method)
##       Live   Video  Both  
## Live         0.0135 <.0001
## Video 0.0135        <.0001
## Both  <.0001 <.0001

Re-test the relationship between body size and individual fighting ability using live-only and video-only data:

# Appropriately subsample 'fishtree':
fishtree_live <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)])
fishtree_live_nosp <- drop.tip(fishtree, c("Spotted Porcupinefish", fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_live)]))
fishtree_video <- drop.tip(fishtree, fishtree$tip.label[!fishtree$tip.label %in% rownames(nojp_scored)])

library(caper)
live_comparative <- comparative.data(data = sheet_phylo_test, phy = fishtree_live, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
live_comparative_nosp <- comparative.data(data = sheet_phylo_test_nosp, phy = fishtree_live_nosp, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)
video_comparative <- comparative.data(data = scored_phylo_test, phy = fishtree_video, names.col = "Species", vcv = TRUE, warn.dropped = TRUE)

live_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
summary(live_pgls_CBI_Mass)
## 
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12134 -0.09154 -0.02659  0.07927  0.14682 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.00051509
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.803, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate  Std. Error t value Pr(>|t|)  
## (Intercept)  1.04492043  0.41860472  2.4962  0.02192 *
## Mass        -0.00048654  0.00074678 -0.6515  0.52251  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09167 on 19 degrees of freedom
## Multiple R-squared: 0.02185, Adjusted R-squared: -0.02963 
## F-statistic: 0.4245 on 1 and 19 DF,  p-value: 0.5225
live_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_CBI_Mass2_3)
## 
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14217 -0.08418 -0.02553  0.07819  0.14423 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.00083993
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.792, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  1.0018370  0.4523589  2.2147   0.0392 *
## Mass2_3     -0.0018632  0.0078794 -0.2365   0.8156  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09256 on 19 degrees of freedom
## Multiple R-squared: 0.002934,    Adjusted R-squared: -0.04954 
## F-statistic: 0.05591 on 1 and 19 DF,  p-value: 0.8156
live_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass)
## 
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.037100 -0.016062  0.006017  0.029263  0.070629 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.663
##    lower bound : 0.000, p = 0.10993
##    upper bound : 1.000, p = 0.01489
##    95.0% CI   : (NA, 0.973)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 0.41376043 0.13109752  3.1561 0.005201 **
## Mass        0.00020690 0.00027338  0.7568 0.458435   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03139 on 19 degrees of freedom
## Multiple R-squared: 0.02926, Adjusted R-squared: -0.02183 
## F-statistic: 0.5728 on 1 and 19 DF,  p-value: 0.4584
live_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3)
## 
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.067679 -0.011879 -0.003656  0.017396  0.048635 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.631
##    lower bound : 0.000, p = 0.13084
##    upper bound : 1.000, p = 0.011061
##    95.0% CI   : (NA, 0.966)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.3593571  0.1397527  2.5714  0.01869 *
## Mass2_3     0.0033572  0.0028881  1.1624  0.25946  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03047 on 19 degrees of freedom
## Multiple R-squared: 0.06639, Adjusted R-squared: 0.01726 
## F-statistic: 1.351 on 1 and 19 DF,  p-value: 0.2595
live_pgls_CBI_Mass_nosp <- pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass_nosp)
## 
## Call:
## pgls(formula = CBI ~ Mass, data = live_comparative_nosp, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.140248 -0.004507  0.026080  0.089520  0.192394 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.00037668
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.787, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.8128104  0.4256586  1.9095  0.07226 .
## Mass        0.0017328  0.0016737  1.0353  0.31424  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0873 on 18 degrees of freedom
## Multiple R-squared: 0.0562,  Adjusted R-squared: 0.003767 
## F-statistic: 1.072 on 1 and 18 DF,  p-value: 0.3142
live_pgls_CBI_Mass2_3_nosp <- pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_CBI_Mass2_3_nosp)
## 
## Call:
## pgls(formula = CBI ~ Mass2_3, data = live_comparative_nosp, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.140775 -0.004609  0.022412  0.089080  0.188876 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 1.000
##    lower bound : 0.000, p = 0.00039991
##    upper bound : 1.000, p = 1    
##    95.0% CI   : (0.785, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.750749   0.450167  1.6677   0.1127
## Mass2_3     0.012323   0.011456  1.0756   0.2963
## 
## Residual standard error: 0.0871 on 18 degrees of freedom
## Multiple R-squared: 0.06039, Adjusted R-squared: 0.008192 
## F-statistic: 1.157 on 1 and 18 DF,  p-value: 0.2963
live_pgls_WinFreq_Mass_nosp <- pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass_nosp)
## 
## Call:
## pgls(formula = WinFreq ~ Mass, data = live_comparative_nosp, 
##     lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.014123  0.002685  0.014639  0.036257  0.064441 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.647
##    lower bound : 0.000, p = 0.079667
##    upper bound : 1.000, p = 0.0081626
##    95.0% CI   : (NA, 0.959)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.31444575 0.14286027  2.2011  0.04102 *
## Mass        0.00116556 0.00070943  1.6430  0.11775  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03031 on 18 degrees of freedom
## Multiple R-squared: 0.1304,  Adjusted R-squared: 0.0821 
## F-statistic: 2.699 on 1 and 18 DF,  p-value: 0.1177
live_pgls_WinFreq_Mass2_3_nosp <- pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp, lambda = "ML"); summary(live_pgls_WinFreq_Mass2_3_nosp)
## 
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = live_comparative_nosp, 
##     lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.061269 -0.013488 -0.002605  0.015575  0.042127 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.643
##    lower bound : 0.000, p = 0.079968
##    upper bound : 1.000, p = 0.0076034
##    95.0% CI   : (NA, 0.957)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 0.2629693  0.1564604  1.6807  0.11008  
## Mass2_3     0.0087316  0.0049639  1.7590  0.09556 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02997 on 18 degrees of freedom
## Multiple R-squared: 0.1467,  Adjusted R-squared: 0.09928 
## F-statistic: 3.094 on 1 and 18 DF,  p-value: 0.09556
video_pgls_CBI_Mass <- pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
summary(video_pgls_CBI_Mass)
## 
## Call:
## pgls(formula = CBI ~ Mass, data = video_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12487 -0.07136 -0.04121  0.01673  0.10692 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.049
##    lower bound : 0.000, p = 0.87521
##    upper bound : 1.000, p = 0.00084113
##    95.0% CI   : (NA, 0.719)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.6174704  0.2932005  2.1060  0.05374 .
## Mass        -0.0014836  0.0020582 -0.7208  0.48287  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07338 on 14 degrees of freedom
## Multiple R-squared: 0.03579, Adjusted R-squared: -0.03308 
## F-statistic: 0.5196 on 1 and 14 DF,  p-value: 0.4829
video_pgls_CBI_Mass2_3 <- pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_CBI_Mass2_3)
## 
## Call:
## pgls(formula = CBI ~ Mass2_3, data = video_comparative, lambda = "ML")
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11810 -0.05457 -0.01404  0.03946  0.10942 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.054
##    lower bound : 0.000, p = 0.86644
##    upper bound : 1.000, p = 0.00091886
##    95.0% CI   : (NA, 0.725)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.6628957  0.3641409  1.8204  0.09013 .
## Mass2_3     -0.0099129  0.0152288 -0.6509  0.52563  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07362 on 14 degrees of freedom
## Multiple R-squared: 0.02938, Adjusted R-squared: -0.03995 
## F-statistic: 0.4237 on 1 and 14 DF,  p-value: 0.5256
video_pgls_WinFreq_Mass <- pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass)
## 
## Call:
## pgls(formula = WinFreq ~ Mass, data = video_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.035535 -0.011008  0.001371  0.023994  0.063956 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.536
##    lower bound : 0.000, p = 0.21789
##    upper bound : 1.000, p = 0.15603
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate  Std. Error t value Pr(>|t|)  
## (Intercept)  0.28375197  0.13948795  2.0342  0.06133 .
## Mass        -0.00043991  0.00083487 -0.5269  0.60650  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02897 on 14 degrees of freedom
## Multiple R-squared: 0.01945, Adjusted R-squared: -0.05059 
## F-statistic: 0.2776 on 1 and 14 DF,  p-value: 0.6065
video_pgls_WinFreq_Mass2_3 <- pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML"); summary(video_pgls_WinFreq_Mass2_3)
## 
## Call:
## pgls(formula = WinFreq ~ Mass2_3, data = video_comparative, lambda = "ML")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.031340 -0.010905 -0.002889  0.022627  0.064360 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.558
##    lower bound : 0.000, p = 0.21225
##    upper bound : 1.000, p = 0.16822
##    95.0% CI   : (NA, NA)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.2890066  0.1631528  1.7714  0.09825 .
## Mass2_3     -0.0025029  0.0059583 -0.4201  0.68081  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02923 on 14 degrees of freedom
## Multiple R-squared: 0.01245, Adjusted R-squared: -0.05809 
## F-statistic: 0.1765 on 1 and 14 DF,  p-value: 0.6808
# Check the distribution of the residuals:
par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass)

par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3)

par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass)

par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3)

par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass_nosp)

par(mfrow=c(2,2))
plot(live_pgls_CBI_Mass2_3_nosp)

par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass_nosp)

par(mfrow=c(2,2))
plot(live_pgls_WinFreq_Mass2_3_nosp)

par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass)

par(mfrow=c(2,2))
plot(video_pgls_CBI_Mass2_3)

par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass)

par(mfrow=c(2,2))
plot(video_pgls_WinFreq_Mass2_3)

Generate dataset-specific versions of Figure 1:

png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_live.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))

par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(live_pgls_CBI_Mass, lwd = 5)
abline(live_pgls_CBI_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass)$r.squared
pvalue <- summary(live_pgls_CBI_Mass)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass_nosp)$coef[2,4]
text(600, 2.95, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 2.75, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 2.5, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 1.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 1.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_CBI_Mass2_3, lwd = 5)
abline(live_pgls_CBI_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(live_pgls_CBI_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_CBI_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_CBI_Mass2_3_nosp)$coef[2,4]
text(58, 2.75, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 2.55, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 2.3, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.6, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.4, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.15, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass, sheet_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass, lwd = 5)
abline(live_pgls_WinFreq_Mass_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass_nosp)$coef[2,4]
text(600, 0.9, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(600, 0.9 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(600, 0.9 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(600, 0.35, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(600, 0.35 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)

par(mar=c(5, 4, 2, 2))
plot(sheet_phylo_test$Mass2_3, sheet_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = ifelse(sheet_phylo_test$Mass == sheet_phylo_test$Mass[sheet_phylo_test$Species == "Spotted Porcupinefish"], "gray70", "black"))
abline(live_pgls_WinFreq_Mass2_3, lwd = 5)
abline(live_pgls_WinFreq_Mass2_3_nosp, lwd = 5, col = "gray70")
lambda <- live_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(live_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(live_pgls_WinFreq_Mass2_3)$coef[2,4]
lambda_nosp <- live_pgls_WinFreq_Mass2_3_nosp$param[2]
pseudor2_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$r.squared
pvalue_nosp <- summary(live_pgls_WinFreq_Mass2_3_nosp)$coef[2,4]
text(58, 0.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(58, 0.5 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(58, 0.5 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
text(58, 0.2, cex = 0.8, bquote(lambda == .(spec_dec(lambda_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2_nosp, 3))), col = "gray70", pos = 4)
text(58, 0.2 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue_nosp, 3))), col = "gray70", pos = 4)
dev.off()

png("/Users/David/Documents/College/2017–18/Winter 2018/FBQ paper/Behav Ecol revisions/Supp_Fig1_video.png", width = 3600, height = 3600, pointsize = 96)
layout(matrix(c(seq(1,4,1)), nrow=2, ncol=2, byrow=TRUE))

par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$CBI, xlab = "Body mass (g)", ylab = "CBI", main = "Isometric scaling", pch = 16, col = "black")
title("(a)", line = -1, adj = 0, outer = TRUE)
title("(b)", line = -1, outer = TRUE)
title("(c)", line = -20, adj = 0, outer = TRUE)
title("(d)", line = -20, outer = TRUE)
abline(video_pgls_CBI_Mass, lwd = 5)
lambda <- video_pgls_CBI_Mass$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass)$r.squared
pvalue <- summary(video_pgls_CBI_Mass)$coef[2,4]
text(150, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)

par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$CBI, xlab = "Body mass raised to 2/3", ylab = "", main = "Allometric scaling", pch = 16, col = "black")
abline(video_pgls_CBI_Mass2_3, lwd = 5)
lambda <- video_pgls_CBI_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_CBI_Mass2_3)$r.squared
pvalue <- summary(video_pgls_CBI_Mass2_3)$coef[2,4]
text(25, 2.5, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 2.3, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 2.05, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)

par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass, scored_phylo_test$WinFreq, xlab = "Body mass (g)", ylab = "Frequency of winning", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass, lwd = 5)
lambda <- video_pgls_WinFreq_Mass$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass)$coef[2,4]
text(150, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(150, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(150, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)

par(mar=c(5, 4, 2, 2))
plot(scored_phylo_test$Mass2_3, scored_phylo_test$WinFreq, xlab = "Body mass raised to 2/3", ylab = "", pch = 16, col = "black")
abline(video_pgls_WinFreq_Mass2_3, lwd = 5)
lambda <- video_pgls_WinFreq_Mass2_3$param[2]
pseudor2 <- summary(video_pgls_WinFreq_Mass2_3)$r.squared
pvalue <- summary(video_pgls_WinFreq_Mass2_3)$coef[2,4]
text(25, 0.55, cex = 0.8, bquote(lambda == .(spec_dec(lambda, 3))), pos = 4)
text(25, 0.55 - 1.2/18, cex = 0.8, bquote(italic(R)^2 == .(spec_dec(pseudor2, 3))), pos = 4)
text(25, 0.55 - 1.2/8, cex = 0.8, bquote(italic(p) == .(spec_dec(pvalue, 3))), pos = 4)
dev.off()

Fit Lanchester’s laws as well as the body size-only and group size-only models to the live-only and video-only datasets:

live_group <- sheet4[sheet4$Group.A != sheet4$Group.B,]
video_group <- scored4[scored4$Group.A != scored4$Group.B,]

# Exclude cornetfish:
live_group <- live_group[live_group$Species.A != "Cornetfish" &
                         live_group$Species.B != "Cornetfish",]

# Number of data points:
nrow(live_group)
## [1] 59
nrow(video_group)
## [1] 24
# Get body size info and use it to test the four models:
lanchester_on_frames <- function(dataframe) {
  n <- ncol(dataframe)
  dataframe$Group.A <- as.numeric(dataframe$Group.A)
  dataframe$Group.B <- as.numeric(dataframe$Group.B)
  for(i in 1:nrow(dataframe)) {
    dataframe[i, n+1] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.A[i])]
    dataframe[i, n+2] <- lengthmass$Mass[lengthmass$Species == as.character(dataframe$Species.B[i])]
    dataframe[i, n+3] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.A[i])]
    dataframe[i, n+4] <- lengthmass$Mass2_3[lengthmass$Species == as.character(dataframe$Species.B[i])]
    if (dataframe$Group.A[i]*dataframe[i, n+1] > dataframe$Group.B[i]*dataframe[i, n+2]) {
      dataframe[i, n+5] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+5] <- dataframe$Species.B[i] 
    }
    if ((dataframe$Group.A[i]^2)*dataframe[i, n+1] > (dataframe$Group.B[i]^2)*dataframe[i, n+2]) {
      dataframe[i, n+6] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+6] <- dataframe$Species.B[i] 
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+5])) == TRUE) {
      dataframe[i, n+7] <- 1
    } else {
      dataframe[i, n+7] <- 0
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+6])) == TRUE) {
      dataframe[i, n+8] <- 1
    } else {
      dataframe[i, n+8] <- 0
    }
    if (dataframe$Group.A[i]*dataframe[i, n+3] > dataframe$Group.B[i]*dataframe[i, n+4]) {
      dataframe[i, n+9] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+9] <- dataframe$Species.B[i] 
    }
    if ((dataframe$Group.A[i]^2)*dataframe[i, n+3] > (dataframe$Group.B[i]^2)*dataframe[i, n+4]) {
      dataframe[i, n+10] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+10] <- dataframe$Species.B[i] 
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+9])) == TRUE) {
      dataframe[i, n+11] <- 1
    } else {
      dataframe[i, n+11] <- 0
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+10])) == TRUE) {
      dataframe[i, n+12] <- 1
    } else {
      dataframe[i, n+12] <- 0
    }
    if (dataframe$Group.A[i] > dataframe$Group.B[i]) {
      dataframe[i, n+13] <- dataframe$Species.A[i]
    } else {
      dataframe[i, n+13] <- dataframe$Species.B[i]
    }
    if (dataframe[i, n+1] > dataframe[i, n+2]) {
      dataframe[i, n+14] <- dataframe$Species.A[i] 
    } else {
      dataframe[i, n+14] <- dataframe$Species.B[i] 
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+13])) == TRUE) {
      dataframe[i, n+15] <- 1
    } else {
      dataframe[i, n+15] <- 0
    }
    if (identical(as.character(dataframe$Winner[i]), as.character(dataframe[i, n+14])) == TRUE) {
      dataframe[i, n+16] <- 1
    } else {
      dataframe[i, n+16] <- 0
    }
  }
  colnames(dataframe)[(n+1):(n+16)] <- c("BodySize.A", "BodySize.B", "BodySize2_3.A",
                                         "BodySize2_3.B", "Linear.Law.Prediction", 
                                         "Square.Law.Prediction", "Linear.True",
                                         "Square.True", "Linear.Law2_3.Prediction",
                                         "Square.Law2_3.Prediction", "Linear2_3.True",
                                         "Square2_3.True", "Group.Size.Only.Prediction",
                                         "Body.Size.Only.Prediction",
                                         "Group.Size.Only.True", "Body.Size.Only.True")
  return(dataframe)
}

live_group <- lanchester_on_frames(live_group)
video_group <- lanchester_on_frames(video_group)

# Proportion of successful predictions:
sum(live_group$Linear.True)/nrow(live_group)
## [1] 0.6949153
sum(live_group$Square.True)/nrow(live_group)
## [1] 0.5254237
sum(live_group$Linear2_3.True)/nrow(live_group)
## [1] 0.6101695
sum(live_group$Square2_3.True)/nrow(live_group)
## [1] 0.4237288
sum(live_group$Group.Size.Only.True)/nrow(live_group)
## [1] 0.2542373
sum(live_group$Body.Size.Only.True)/nrow(live_group)
## [1] 0.6779661
sum(video_group$Linear.True)/nrow(video_group)
## [1] 0.375
sum(video_group$Square.True)/nrow(video_group)
## [1] 0.08333333
sum(video_group$Linear2_3.True)/nrow(video_group)
## [1] 0.1666667
sum(video_group$Square2_3.True)/nrow(video_group)
## [1] 0.04166667
sum(video_group$Group.Size.Only.True)/nrow(video_group)
## [1] 0
sum(video_group$Body.Size.Only.True)/nrow(video_group)
## [1] 0.75
# Test significance:
dbinom(sum(live_group$Linear.True), size = nrow(live_group), prob = 1/2) 
## [1] 0.001123269
dbinom(sum(live_group$Square.True), size = nrow(live_group), prob = 1/2)
## [1] 0.09596023
dbinom(sum(live_group$Linear2_3.True), size = nrow(live_group), prob = 1/2) 
## [1] 0.02501637
dbinom(sum(live_group$Square2_3.True), size = nrow(live_group), prob = 1/2)
## [1] 0.05253438
dbinom(sum(live_group$Group.Size.Only.True), size = nrow(live_group), prob = 1/2) 
## [1] 6.920778e-05
dbinom(sum(live_group$Body.Size.Only.True), size = nrow(live_group), prob = 1/2)
## [1] 0.002423897
dbinom(sum(video_group$Linear.True), size = nrow(video_group), prob = 1/2) 
## [1] 0.07793331
dbinom(sum(video_group$Square.True), size = nrow(video_group), prob = 1/2)
## [1] 1.645088e-05
dbinom(sum(video_group$Linear2_3.True), size = nrow(video_group), prob = 1/2) 
## [1] 0.000633359
dbinom(sum(video_group$Square2_3.True), size = nrow(video_group), prob = 1/2)
## [1] 1.430511e-06
dbinom(sum(video_group$Group.Size.Only.True), size = nrow(video_group), prob = 1/2) 
## [1] 5.960464e-08
dbinom(sum(video_group$Body.Size.Only.True), size = nrow(video_group), prob = 1/2)
## [1] 0.008022547

How many interactions involving damselfish do we have in the live-only dataset and how many in the video-only dataset?

nrow(live_group[live_group$Species.A == "Damselfish" | live_group$Species.B == "Damselfish",])/nrow(live_group)
## [1] 0.3559322
nrow(video_group[video_group$Species.A == "Damselfish" | video_group$Species.B == "Damselfish",])/nrow(video_group)
## [1] 0.7916667

References